You are right:.

DOp[v][f0x[t_, x_, v_] = 
   E^(-(1 + 1/2) v t - x^2/2) HermiteH[1, x] + 
    2/10 E^(-(4 + 1/2) v t - x^2/2) HermiteH[4, x] + 
    1/100 E^(-(6 + 1/2) v t - x^2/2) HermiteH[6, x]] // FullSimplify

After making v a parameter to  DOp  and defining the start distribution f0x,  forgot to  update the definition

  nds[p_, {t_, x_}, v_, {t_, 0, T_}, {x_, a_, b_}] :=
 NDSolve[
  {DOp[v][p[t, x]] == 0,
   p[0, x] == f0x[0, x, v],
   Derivative[0, 1][p][t, a] == 0, Derivative[0, 1][p][t, b] == 0},
  p[t, x], {t, 0, T}, {x, a, b}]

nds[p, {t, x}, 1, {t, 0, 1}, {x, -4, 4}]


q[t_, x_] = 
 p[t, x] /. {nds[p, {t, x}, 1, {t, 0, 1}, {x, -4, 4}][[1]]} // Quiet


Manipulate[
  Plot[{q[t, x], f0x[t, x, 1]}, {x, -4, 4}, 
    PlotStyle -> {{Red, Thickness[0.02]}, {Black}}, 
    PlotRange -> {-6, 6}], {t, 0, 1}]

Hi Athanasios, seems I dont quite understand..

LegrendreP[n,x] is orthogonal on (-1,1) with measure (1-x^2), HermiteH[n,x ] is orthogonal on -inf, inf with measure e^(-x^2).
e^(-x^2) H_n(x) is orthogonal on -inf, inf with measure 1 e^(-x^2) Hn(x) are eigenfunctions of -dxx +x^2 with eigenvalue 2n+1

It follows without any approximation that

Sum[a_n  e^(-(2n+1)t  HermiteH[n,x]) ,{ n, 0, inf}]

solves

d_t == d_xx +  x^2 

Yor start distribution is of that form for some a_n

Of course you may approximate e^-x2/2 HermiteH[n,x] by LegendreP[m,x] on (-1.1), But only using the measure (1-x^2) dx and normailzation to an orthonormal basis will yield the correect approximation.

Regards Roland

In[1]:= 
DOp[v_] = ((D[#, t] - v/2 D[#, {x, 2}] + v/2  x^2 #)  &);

In[3]:= Assuming[n \[Element] PositiveIntegers , 
 DOp[v][E^(-(n + 1/2) v t - x^2/2)  HermiteH[n, x]] // FullSimplify]

Out[3]= 0

In[5]:=  
DOp[v][f0x[t_, x_, v_] = 
   E^(-(1 + 1/2) v t - x^2/2)  HermiteH[1, x] + 
    2/10 E^(-(4 + 1/2) v t - x^2/2)  HermiteH[4, x] +
    1/100  E^(-(6 + 1/2) v t - x^2/2)  HermiteH[6, x] ] // FullSimplify

Out[5]= 0

In[6]:= nds[p_, {t_, x_}, v_, {t_, 0, T_}, {x_, a_, b_}] :=
 NDSolve[   {DOp[ p[t,x] ==0,   
   p[0, x] ==  2 E^(-(x^2/2)) x + 1/5 E^(-(x^2/2)) (12 - 48 x^2 + 16 x^4) + 
     1/100 E^(-(x^2/2)) (-120 + 720 x^2 - 480 x^4 + 64 x^6)       
   Derivative[0, 1][p][t, a] == 0, Derivative[0, 1][p][t, b] == 0}, 
  p[t, x], {t, 0, T}, {x, a, b}]

In[7]:= 
q[t_, x_] = p[t, x] /. {nds[p, {t, x}, 1, {t, 0, 1}, {x, -4, 4}][[1]]}

During evaluation of In[7]:= NDSolve::ibcinc: Warning: boundary and initial conditions are inconsistent.

Out[7]= {InterpolatingFunction[{{0., 1.}, {-4., 
   4.}}, {5, 5, 1, {76, 75}, {4, 6}, 0, 0, 0, 0, Automatic, {}, {}, 
    False}, {CompressedData[s=
"]}, {Automatic, Automatic}][t, x]}

Manipulate[
 Plot[{q[t, x], f0x[t, x, 1]}, {x, -4, 4}, 
  PlotStyle -> {{Red, Thickness[0.02]}, {Black}}, 
  PlotRange -> {-6, 6}], {t, 0, 1}]

Just a question: You are using the LegendreP expansion together with the Gaussian measure, but the eigenfunctions on the infinite interval are Parabolic cylinder functions with Hermite polynomials as the special case of normalizable functions on all of R.

The LegendreP is of no use outside (-1,1), What is the intention of their use here?

Regards Roland