(* Source code written in Mathematica 6.0 by Steve Byrnes, Feb. 2011. This source code is public domain. *)
(* Shows classical and quantum trajectory animations for a harmonic potential. Assume m=w=hbar=1. *)
ClearAll["Global`*"]
(*** Wavefunctions of the energy eigenstates ***)
psi[n_, x_] := (2^n*n!)^(-1/2)*Pi^(-1/4)*Exp[-x^2/2]*HermiteH[n, x];
energy[n_] := n + 1/2;
psit[n_, x_, t_] := psi[n, x] Exp[-I*energy[n]*t];
(*** A random time-dependent state ***)
SeedRandom[1];
CoefList = Table[Random[]*Exp[2 Pi I Random[]], {n, 0, 4}];
CoefList = CoefList/Norm[CoefList];
Randpsi[x_, t_] := Sum[CoefList[[n + 1]]*psit[n, x, t], {n, 0, 4}];
(*** A coherent state (or "Glauber state") ***)
CoherentState[b_, x_, t_] := Exp[-Abs[b]^2/2] Sum[b^n*(n!)^(-1/2)*psit[n, x, t], {n, 0, 15}];
(*** Make the classical plots...a red ball anchored to the origin by a gray spring. ***)
classical1[t_, max_] := ListPlot[{{max Cos[t], 0}}, PlotStyle -> Directive[Red, AbsolutePointSize[15]]];
zigzag[x_] := Abs[(x + 0.25) - Round[x + 0.25]] - .25;
spring[x_, left_, right_] := (.9 zigzag[3 (x - left)/(right - left)])/(1 + Abs[right - left]);
classical2[t_, max_] := Plot[spring[x, -5, max Cos[t]], {x, -5, max Cos[t]}, PlotStyle -> Directive[Gray, Thick]];
classical3 = ListPlot[{{-5, 0}}, PlotStyle -> Directive[Black, AbsolutePointSize[7]]];
classical[t_, max_, label_] := Show[classical2[t, max], classical1[t, max], classical3,
   PlotRange -> {{-5, 5}, {-1, 1}}, Ticks -> None, Axes -> {False, True}, PlotLabel -> label, AxesOrigin -> {0, 0}];
(*** Put all the plots together ***)
SetOptions[Plot, {PlotRange -> {-1, 1}, Ticks -> None, PlotStyle -> {Directive[Thick, Blue], Directive[Thick, Pink]}}];
MakeFrame[t_] := GraphicsGrid[
   {{classical[t + 2, 1.5, "A"], classical[t, 3, "B"]},
    {Plot[{Re[psit[0, x, t]], Im[psit[0, x, t]]}, {x, -5, 5}, PlotLabel -> "C"],
     Plot[{Re[psit[1, x, t]], Im[psit[1, x, t]]}, {x, -5, 5}, PlotLabel -> "D"]},
    {Plot[{Re[psit[2, x, t]], Im[psit[2, x, t]]}, {x, -5, 5}, PlotLabel -> "E"],
     Plot[{Re[psit[3, x, t]], Im[psit[3, x, t]]}, {x, -5, 5}, PlotLabel -> "F"]},
    {Plot[{Re[Randpsi[x, t]], Im[Randpsi[x, t]]}, {x, -5, 5}, PlotLabel -> "G"],
     Plot[{Re[CoherentState[1, x, t]], Im[CoherentState[1, x, t]]}, {x, -5, 5}, PlotLabel -> "H"]}
    }, Frame -> All, ImageSize -> 300];
output = Table[MakeFrame[t], {t, 0, 4 Pi*96/97, 4 Pi/97}];
SetDirectory["C:\Users\Steve\Desktop"]
Export["test.gif", output]

psit[n_, x_, t_] := psi[n, x] Exp[-I*energy[n]*t];