Series expansion in terms of Hermite polynomials

For polynomials, you don't need to do any integrals to find the expansion. Take a polynomial p and a list basis containing the basis functions. Then define a function that takes these two, identifies the variable x, and solves for the coefficients in basis that make the two polynomials equal in terms of their CoefficientLists:

expandPoly[p_, basis_, x_] :=
# /. First@Solve[CoefficientList[#.basis, x] == #2, #] &[
     Array["a", Length[#]], #] & @
     PadRight[CoefficientList[p, x], Length[basis]]

expandPoly[1 + x + 3 x^2 + 7 x^3, HermiteH[Range[4] - 1, x], x]

(* ==> {5/2, 23/4, 3/4, 7/8} *)

Edit

In response to belisarius: if you already know that you're only interested in a basis of HermiteH, you could incorporate that into the function and do away with the specification of the variable basis as follows:

expandPoly[p_, x_] := # /. 
      First @ Solve[
        CoefficientList[#.HermiteH[Range[Length[#]] - 1, x], 
          x] == #2, #] &[Array["a", Length[#]], #] & @ CoefficientList[p, x]

expandPoly[1 + x + 3 x^2 + 7 x^3, x]

(* ==> {5/2, 23/4, 3/4, 7/8} *)

Edit 2

With the general function given as the first solution above, you can specify any set of polynomials that is known to form a basis for degree n or larger. This means the basis functions don't have to be orthogonal polynomials at all.

The Hermite polynomials are orthogonal with respect to the inner product $$\langle f,g \rangle = \int_{-\infty}^{\infty} f(x)g(x)e^{-x^2} \, \mathrm dx.$$ Thus, the $n$-th coefficient can be computed using the inner product of your polynomial with the $n$-th normalized Hermite polynomial.

Example:

p[x_] = 1 + x + x^2 + x^3;
coeffs = Table[
  Integrate[HermiteH[n, x]*p[x]*Exp[-x^2], {x, -∞, ∞}]/
   Integrate[HermiteH[n, x]^2*Exp[-x^2], {x, -∞, ∞}],
  {n, 0, 3}]
(* Out: {3/2, 5/4, 1/4, 1/8} *)

coeffs.Table[HermiteH[n, x], {n, 0, 3}] // Expand
(* Out: 1 + x + x^2 + x^3 *)

The inner product for the Hermite polynomials, $$\langle f, g\rangle \int_{-\infty}^{\infty} f(x)\,g(x)\,e^{-x^2}\;dx\,,$$ has nice formulas for power functions (where $n=a+b$) and for the Hermite polynomials: $$ \begin{align} \langle x^a, x^b \rangle = \langle x^n, 1\rangle &= \frac{1}{2} \left((-1)^n+1\right)\, \Gamma \left(\frac{n+1}{2}\right)\cr \langle H_n(x), H_n(x) \rangle &= \sqrt{\pi}\,2^n n! \cr \end{align} $$

These can be used to give a quick change of basis function for polynomials.

hermiteIP[f_, g_, x_] := With[{coeff = CoefficientList[f g, x]},
   coeff.Table[1/2 (1 + (-1)^(-1 + n)) Gamma[n/2], {n, Length@coeff}]];

hermiteExpand[poly_, var_] /; PolynomialQ[poly, var] :=
 Sum[hermiteIP[poly, HermiteH[n, var], var] H[n, var]/(Sqrt[Pi] 2^n n!),
  {n, 0, Exponent[poly, var]}]

I used H[n, x] as a place holder for HermiteH[n, x].

hermiteExpand[(1 + x)^5, x]
(* 39/4 H[0, x] + 95/8 H[1, x] + 25/4 H[2, x] + 15/8 H[3, x] + 
    5/16 H[4, x] + 1/32 H[5, x]  *)

hermiteExpand[(1 + x)^5, x] /. H -> HermiteH
(* 39/4 H[0, x] + 95/8 H[1, x] + 25/4 H[2, x] + 15/8 H[3, x] + 
    5/16 H[4, x] + 1/32 H[5, x]  *)

% // Factor
(* (1 + x)^5 *)