Zeros of high degree polynomials

Don't use machine numbers, as subtractive cancellation will cause enormous precision loss, as is common with high order polynomials.

You can either work with exact results using Solve:

HermiteH[18, x /. Solve[HermiteH[18,x]==0,x,Reals]] //Simplify

{0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0}

Or you can use the WorkingPrecision option:

HermiteH[18, x /. NSolve[HermiteH[18,x]==0, x, Reals, WorkingPrecision->50]]

{0.*10^-33, 0.*10^-35, 0.*10^-36, 0.*10^-37, 0.*10^-38, 0.*10^-38, 0.*10^-39, 0.*10^-39, 0.*10^-40, 0.*10^-40, 0.*10^-39, 0.*10^-39, 0.*10^-38, 0.*10^-38, 0.*10^-37, 0.*10^-36, 0.*10^-35, 0.*10^-33}

You can see that subtractive cancellation causes a precision loss of about 17 digits for the first and last roots.

For polynomials like HermiteH, these roots are represented in Mathematica as infinite-precision Root objects. The $k$th root of the $n$th Hermite polynomial is, with infinite precision, represented by

R[n_, k_] := Root[HermiteH[n, #] &, k]

What Carl's use of Solve does is simply to make a list of such Root objects. You can work with these objects analytically (using RootReduce etc.), or you can convert them to numerical.

For example, the 7th root of $H_{18}$ would be

r = R[18, 7]
(*    a root around -1.30...    *)

Numerically:

N[r]
(*    -1.30092    *)
N[r, 100]
(*    -1.300920858389617365666265554392610580218134639661226522772309775882782630084141194539623631652544514    *)

analytic transformations:

r^2 // RootReduce
(*    a root around 1.69...    *)
HermiteH[18, r] // RootReduce
(*    0    *)

At least for this case, one can also consider getting the eigenvalues of the Jacobi matrix associated with the Hermite polynomials. Recall that these matrices are constructed from the coefficients of the three-term recurrences that generate the corresponding orthogonal polynomial. Applied to this case, we have:

With[{n = 18}, 
     s1 = Sort[Eigenvalues[N[SparseArray[{{j_, k_} /; Abs[j - k] == 1 :> Sqrt[Min[j, k]/2]},
                                         {n, n}]]]]]
   {-5.04836, -4.24812, -3.57377, -2.96138, -2.3863, -1.83553, -1.30092, -0.776683,
    -0.258268, 0.258268, 0.776683, 1.30092, 1.83553, 2.3863, 2.96138, 3.57377, 4.24812,
    5.04836}

For comparison purposes, let's compare that and a machine-precision evaluation of NSolve[] with an arbitrary-precision evaluation:

s2 = Sort[x /. NSolve[HermiteH[18, x], x]];

sb = Sort[x /. NSolve[HermiteH[18, x], x, WorkingPrecision -> 30]];

{Norm[sb - s1, ∞], Norm[sb - s2, ∞]}
   {4.44089*10^-15, 2.40696*10^-13}

and we see that the Jacobi-based method gives better results.