(*Cell 1. We are interested in the general relationships here; set constants to 1. Run this cell.*)
m:=1
ω:=1
ℏ:=1
b:=(Sqrt[(ℏ/(m ω))])

(*Cell 2. These equations are already written for you. Review what is written here, because you'll use it later. Run this cell.*)
U[x_]:=(1/2)m ω^2 x^2
ψ1[x_]:=A1 Exp[-x^2/(2b^2)]
ψ2[x_]:= A2(x/b)Exp[-x^2/(2b^2)]
ψ3[x_]:=A3 (1-((2x^2)/(b^2)))Exp[-x^2/(2b^2)]
En[n_]:=(n-(1/2))ℏ ω

(*Cell 3. Input the expression for A1 that you found by hand. The other amplitudes are already written for you, in terms of A1*)
A1:=
A2:=(Sqrt[2])*A1
A3:=-A1/(Sqrt[2]) (*The minus sign is to match the wavefunction in book; remember wavefunction is not unique, would create same probability distribution either way.*)

(*Cell 4. Now use the "Integrate" function to check that your normalization constant is correct. The general command Integrate[f,{x,xmin,xmax}] gives the integral of the function f(x) from xmin to xmax. You may type the word "Infinity" into Mathematica and it will understand. You can find the necessary Greek symbols yourself or just copy and paste from above if needed. *)

(*Cell 5. Run this cell to view a plot of the potential and the first 3 wavefunctions.*)
Plot[{U[x],ψ1[x],ψ2[x],ψ3[x]},{x,-3b,3b},PlotRange-> {-0.8,3.2},PlotStyle-> {Dashed, Solid, Solid, Solid},AxesLabel->{"x","U(x) dashed, ψ(x) solid"}]

(*Cell 6. It's pretty, but it's kind of hard to read. How can you *modify* the previous command to plot the functions the way they appear in your book, spaced out vertically according to energy level? Try it here.*)
Plot[{U[x],ψ1[x],ψ2[x],ψ3[x]},{x,-3b,3b},PlotRange-> {-0.8,3.2},PlotStyle-> {Dashed, Solid, Solid, Solid},AxesLabel->{"x","U(x) dashed, ψ(x) solid"}]

(*Cell 7. The general solution to the Schrodinger equation for the quantum harmonic oscillator involves Hermite polynomials. Run the following cell to allow you to view the shape of the probability density for various n. This code is part of the Wolfram Demonstrations Project,contributed by Michael Trott.*)
ψgen[n_,x_]:=1/Sqrt[Sqrt[Pi] 2^n n!] Exp[-x^2/2] HermiteH[n,x]
Manipulate[
Plot[Evaluate[ψgen[n,x]^2],Evaluate@{x,-Sqrt[n]-4,Sqrt[n]+4},Filling->Axis,Frame->True,Axes->False,ImageSize->{500,400},ImagePadding->{{30,20},{20,20}}],
{{n,0,"energy level"},0,30,1},SaveDefinitions->True]

(*Cell 8. Let's see the time-dependence of these wavefunctions. Write a formula for the time-dependent part of the wavefunction, as a function of n, quantum number. The other equations are written for you. To get the imaginary number "i", You can type Esc+ii+Esc and the result will look like this: I. You can also copy-and-paste that I I just gave you.*)
ϕ[n_,t_]:=
Ψ1[x_,t_]:=ψ1[x]*ϕ[1,t]
Ψ2[x_,t_]:=ψ2[x]*ϕ[2,t]
Ψ3[x_,t_]:=ψ3[x]*ϕ[3,t]

(*Cell 9. Run this cell to animate the real (solid) and imaginary (dashed) parts of the wave function*)
Animate[Show[Plot[{Re[Ψ1[x,t]],Im[Ψ1[x,t]],Abs[Ψ1[x,t]]},{x,-3b,3b},PlotRange-> {-2,2},PlotStyle->{Dashed,Dashed,Solid},AxesLabel->{"x","Ψ1"}]],{t,0,4Pi}]
Animate[Show[Plot[{Re[Ψ2[x,t]],Im[Ψ2[x,t]],Abs[Ψ2[x,t]]},{x,-3b,3b},PlotRange-> {-2,2},PlotStyle->{Dashed,Dashed,Solid},AxesLabel->{"x","Ψ2"}]],{t,0,4Pi}]
Animate[Show[Plot[{Re[Ψ3[x,t]],Im[Ψ3[x,t]],Abs[Ψ3[x,t]]},{x,-3b,3b},PlotRange-> {-2,2},PlotStyle->{Dashed,Dashed,Solid},AxesLabel->{"x","Ψ3"}]],{t,0,4Pi}]