U[x_] := x^2/2

p[x_]:=Sqrt[
        2*(e - U[x])]

turn=Solve[U[x] == e, 
        x]

xleft = x /. turn[[1]]

xright = x /. turn[[2]]

level[n_] := 
      Graphics[Line[{{xleft, e}, {xright, e}} /. {e -> 
                n + 1/2}]]

uplot := Plot[U[x],  {x, -3, 3}]

uplot

uplot := 
      Plot[U[x],  {x, -3.1, 3.1}, PlotStyle -> RGBColor[0, 0, 1],  
        DisplayFunction -> Identity]

Show[uplot, 
      Table[level[n],  {n, 0, 
          4}],   DisplayFunction -> $DisplayFunction]

(*now we get the wavefunctons of the HO*) 
?Hermit*

psi[n_Integer, z_] := 
      1/Pi^(1/4) 1/Sqrt[2^n*n!]*Exp[-z^2/2]*
        HermiteH[n, z]

(∫_(-∞)^∞)(( psi(1, 
              z))^2) \[DifferentialD]z

(∫_(-∞)^∞) psi[
              Random[Integer,  {0, 7}], z]^2 \[DifferentialD]z

ri := Random[Integer,  {0, 7}]

Integrate[
      psi[ri, z]*psi[ri, z],  {z, -Infinity,  Infinity}]

Clear[p]; 
    PsiCoh[z_, t_] := 
      1/Pi^(1/4)*
Exp[I p[t]*z - 1/2*(z - x[t])^2]*
        Exp[-I t/2 - 1/2*I*p[t]*x[t]]

x[t_] := 
      Sqrt[2  e]*Cos[t]

    p[t_] := -Sqrt[2  e]*Sin[t]

Clear[e, t]; 
    Plot[Abs[PsiCoh[z, t]] /. {e -> 8,  t -> 0},  {z,  -8, 
        8},  PlotRange -> {0, 1/Pi^(1/4)}]

red = RGBColor[1, 0, 0]; green = RGBColor[0, 1, 0];  
    blue = RGBColor[0, 0, 1];

Plot[{Abs[PsiCoh[z, t]]^2 /. {e -> 8,  
            t -> 0}, 
Abs[
              PsiCoh[z, t]]^2 /. {e -> 8,  t -> Pi/2}, 
        Abs[PsiCoh[z, t]]^2 /. {e -> 8,  
            t -> Pi}},  {z,  -8, 8},  
      PlotRange -> {0, 1/Pi^(1/2)},  
      PlotStyle -> {red, green, blue}]

Integrate[Abs[PsiCoh[z, t]]^2,  {z, -Infinity, Infinity},  
      Assumptions -> {t > 0, e > 0}]

a[n_] := 
      FullSimplify[
        Exp[I *(n + 1/2)*t] 
          Integrate[psi[n, z]*PsiCoh[z, t],  {z, -Infinity, Infinity},  
            Assumptions -> {n > 0,  t > 0,  e > 0,  Cos[t] > 0,  
                Sin[t] > 0}],  e > 0]

FullSimplify[a[1] /. t -> 0,  e > 0]

Abs[a[1]]^2 /. {t -> Random[], e -> 4. }

pois[n_Integer, nav_] := E^(-nav)*nav^n/n!

NSum[pois[n, nav],  {n, 0, Infinity}] /. nav -> ri

pois[1, 4. ]

aa[n_] := 
      FullSimplify[
        Integrate[psi[n, z]*PsiCoh[z, 0],  {z, -Infinity, Infinity},  
          Assumptions -> {n >= 0,  e > 0}],  e > 0]

aa[1]

e = 8. ; Timing[list = Table[aa[n],  {n, 0, 5}]]

PsiApp[z_, t_, k_] := 
      Sum[aa[n]*Exp[-i (n + 1/2)*t]*psi[n, z],  {n, 0, 
          k}]

ploApp[T_, kk_] := 
      Block[{t = T, k = kk},  
        Plot[Evaluate[Abs[PsiApp[z, t, k]]^2],  {z, -8, 8}, 
          PlotRange -> {0,  .7},  
          PlotStyle -> {red, Dashing[{ .2/kk}]}, 
          DisplayFunction -> Identity]]

ploCoh[T_] := 
      Block[{t = T},  
        Plot[Evaluate[Abs[PsiCoh[z, t]]^2],  {z, -8, 8}, 
          PlotRange -> {0,  .7},  PlotStyle -> blue, 
          DisplayFunction -> Identity]]

Show[ploApp[0, 4], ploApp[0, 8], ploApp[0, 16],  ploCoh[0],  
      DisplayFunction -> $DisplayFunction]