(* Disclaimer: zero efforts have been made to optimize anything in the notebook below *)
(*cutoff parameter *) nn=100;
(* dimensionless interaction strength *) μ=0.1;(* Hamiltonian matrix *)ham=Table[(Sqrt[n(n-1)(n-2)(n-3)]*KroneckerDelta[m,n-4]+(4n-2)*Sqrt[n(n-1)]*KroneckerDelta[m,n-2]+(6n^2+6n+3)*KroneckerDelta[m,n]+Sqrt[(n+1)(n+2)]*(4n+6)*KroneckerDelta[m,n+2]+Sqrt[(n+1)(n+2)(n+3)(n+4)]*KroneckerDelta[m,n+4])*μ+(n+1/2)*KroneckerDelta[m,n],{n,0,nn},{m,0,nn}];
(* eigenvalues and eigenvectors *)res=Eigensystem[ham];
evals=res[[1]];
evecs=res[[2]];

(* energy eigenvalues *)Table[evals[[j]],{j,1,nn+1}]

(* overlaps <(Subscript[ψ, j])|0> *) Table[Chop[evecs[[j,1]]],{j,1,nn+1}]

(* Harmonic oscillator wave functions in dimensionless variables *)
psi[n_,x_]:=1/Sqrt[Factorial[n]*2^n]*HermiteH[n,x/Sqrt[2]]*Exp[-x^2/4]/(2*Pi)^(1/4)
(* dimensionless probability density *)p[z_,tau_]:=Abs[Sum[psi[k-1,z]*Exp[-I*tau*evals[[n]]]*evecs[[n,1]]*evecs[[n,k]],{k,1,nn+1,2},{n,1,nn+1,2}]]^2
(* here we sum only over half the modes because the other overlaps are zero due to parity symmetry *)

(* test that we recover the t=0 result *)Show[{ListPlot[Table[{0.1*j,p[0.1*j,0]},{j,-40,40}]],Plot[Abs[psi[0,x]]^2,{x,-5,5}]}]

res05=Table[{0.1*j,p[0.1*j,0.5]},{j,-30,30}];
res1=Table[{0.1*j,p[0.1*j,1.0]},{j,-30,30}];
res15=Table[{0.1*j,p[0.1*j,1.5]},{j,-30,30}];
res2=Table[{0.1*j,p[0.1*j,2]},{j,-30,30}];

ListPlot[Table[{0.1*j,p[0.1*j,1.0]},{j,-30,30}],PlotRange->All]

ListPlot[{res05,res1,res15,res2}]