**First point** : in Maple the Hermite polynomials are mutually orthogonal with respect to the measure **exp(-y^2) dy**; (which differs from their definition in Statistics where there are mutually orthogonal with respect to the measure** df(y),** where **f** denotes the PDF of a standard gaussian random variable, see for instance Hermite_polynomials).

~~Accordingly your definition of ~~**HN** is wrong, and even doubly wrong because of the unjustified presence of an extra **sqrt**.

You should have define **HN** this way

~~HN := (n, y) -> exp(-y^2)*simplify(HermiteH(n, y)) / (sqrt(Pi)*2^n*n!):
~~

**Second point**: in a reply @acer suggested you to use **[fsolve(simplify(HermiteH(7, x)))]; ** instead of ** evalf(allvalues(RootOf(HermiteH(7, x), x)));** and you answered him "*One solution is not enough, that's why I tried with solve. I need to get all the possible solutions, then discard complex ones and study the behavior of the ones that only provide pure real coefficients for the expansion".*

**It's pure nonsense! **

All the zeros of Hermite polynomials are real. For a quick and simple proof look here the-roots-of-hermite-polynomials-are-all-real

What you want to compute are the values named "Gauss points", or "Gauss integration points", which are indeed the zeros of the Hermite polynomials..Have you ever heard of a Gauss quadrature method where the integration points would be complex ?

These zeros have symmetric values, implying that Hermite polynomials of odd degrees have 0 as a root.

You should follow acer's advice, all the more that the way you compute the Gauss points gives wrong results:

r := evalf(allvalues(RootOf(HermiteH(7, x), x)));
0.8162878829, -1.673551629, -0.8162878829, 0.

**Only 4 roots instead of 7!**

(which explains why @tomleslie's solution contains so many undeterminate coefficients)

Here are the "corrrect" values

fsolve( simplify(HermiteH(7, x), 'HermiteH') );
-2.651961357, -1.673551629, -0.8162878829, 0., 0.8162878829, 1.673551629, 2.651961357

Another way to get them is

evalf([solve(simplify(HermiteH(7, x), 'HermiteH'))]);
sort(Re~(%));
[ -11
[0., 2.651961356 + 5.356554700 10 I,
-11
-2.651961356 - 5.356554700 10 I,
-10
0.8162878825 + 2.479096188 10 I,
-10
-0.8162878825 - 2.479096188 10 I,
-10
1.673551629 - 3.850404915 10 I,
-10 ]
-1.673551629 + 3.850404915 10 I]
[-2.651961356, -1.673551629, -0.8162878825, 0., 0.8162878825, 1.673551629, 2.651961356]

**Finally**: Here is a full computation of the *ansatz* **u7**

Download HermiteH_.mw

Oldies but Goodies: jresv48n2p111_A1b.pdf