Time Dependent Schrödinger Equation
T J Atherton, Tufts University

How to use this notebook:
1) Run the whole notebook by choosing Edit>Select All and then Shift-Enter.
2) Choose a potential in the “Prepare a potential” section of the notebook.
3) Construct a wavepacket in the “Prepare a wavepacket” section of the notebook.
4) Go to the interactive video at the bottom and click “Go” to start the video. 
5) Modify the potential or wavepacket as you like and click “Update” to refresh the display.

Select region of space and discretization

Region of space to solve over

xsolveDomain=N[{-40,40}];

Region of space to visualize solution over

domainWidth=N[(xsolveDomain[[2]]-xsolveDomain[[1]])];

domainCenter = (xsolveDomain[[2]]+xsolveDomain[[1]])/2;

xvisDomain=N[{domainCenter-domainWidth/8,domainCenter+domainWidth/8}]

Discretization choices

δx=N[domainWidth/1000];

δt = N[δx];

tmax= 100000 δt;

pause = 0.02;

Prepare solver

(*PPinfwell=-1/2 SparseArray[Join[Table[{i,i}-> -2,{i,NN}],Table[{Wrap[i-1],i}-> 1,{i,2,NN}],Table[{Wrap[i+1],i}-> 1,{i,1,NN-1}]]]/δx^2;*)

joinCoordinates[xx_,f_]:= Transpose[{xx,f}];

Wrap[i_,NN_]:=Which[i<1,NN,i>NN,1,True,i];

prepareSolver[V_,ψ0_,domain_,δx_,δt_]:= Block[{xx,NN,vv,ps0,VVx,II,PPperiodic,HH,h1,h2,Lh2},
xx=Table[x,{x,domain[[1]],domain[[2]],δx}];
NN=Length[xx];
     vv=Table[V[x],{x,domain[[1]],domain[[2]],δx}];
ps0=Table[ψ0[x,x0,k0,σk],{x,domain[[1]],domain[[2]],δx}];
(*ps0=ps0/Sqrt[ps0.Conjugate[ps0]];*)
    VVx=SparseArray[Table[{i,i}-> vv[[i]],{i,NN}]];
    II=SparseArray[Table[{i,i}-> 1,{i,NN}]];
    
PPperiodic=-1/2 SparseArray[Join[Table[{i,i}-> -2,{i,NN}],Table[{Wrap[i-1,NN],i}-> 1,{i,NN}],Table[{Wrap[i+1,NN],i}-> 1,{i,NN}]]]/δx^2;
HH=PPperiodic+VVx;
    h1=(II- I HH/2 δt);
    h2=(II+ I HH/2 δt);
    Lh2=LinearSolve[h2];
Return[{h1,Lh2,xx,ps0,vv}];
]

Prepare a potential

Use the GUI below to choose a potential function A V(B(x-(Subscript[x, 0]))); you can adjust the constants A, B and (Subscript[x, 0]) to scale and translate the function as you like. You can also define an arbitrary potential by modifying these functions—the (Subscript[V, custom]) potential is provided specifically for this purpose.

Vfree[x_]:=0

Vho[x_]:=x^2

Vstep[x_]:= If[x<0,0,1]

Vbarrier[x_]:=Which[0<x<1,1,True,0];

Vcrystal[x_]:= -Which[x>0 && EvenQ[IntegerPart[ x]],1,True,0];

Vcustom[x_]:= 0

Manipulate[V=A VV[B(#-xv0)]&;Plot[A VV[B(x-xv0)],{x,xvisDomain[[1]],xvisDomain[[2]]},PlotRange->{-10,10},PlotStyle-> {Thick,Black},PlotPoints-> 200],{{VV,Vstep},{Vfree,Vstep,Vbarrier,Vcrystal,Vho,Vcustom}},{{A,1},-10,10},{{B,1},0,2},{{xv0,0},xvisDomain[[1]],xvisDomain[[2]]},LocalizeVariables->False,SaveDefinitions-> True]

Prepare a wavepacket

Prepare a starting wavefunction — a normalized Gaussian wavepacket with center x0 and momentum k0 with width in momentum space k0

ψ0[x_,x0_,k0_,σ_]:=(E^((1/4) (x-x0) (4 I k0+(-x+x0) (σ^2)))) (Sqrt[π]) σ/Sqrt[(Sqrt[2])  (π^(3/2)) σ];

Manipulate[Plot[{Re[ψ0[x,x0,k0,σk]],Im[ψ0[x,x0,k0,σk]],Abs[ψ0[x,x0,k0,σk]]},{x,xsolveRange[[1]],xsolveRange[[2]]},PlotRange-> {-1,1},PlotStyle-> {Cyan,Pink,Black}],{{x0,0},xvisDomain[[1]],xvisDomain[[2]]},{{k0,0},0,5},{{σk,0.05},1*^-4,5},LocalizeVariables-> False,SaveDefinitions-> True]

See the Movie!

{h1,Lh2,xx,psi0,vv}=prepareSolver[V,ψ0,xsolveDomain,δx,δt];psi2=psi0;n=0;go=False;Grid[{{Grid[{{Button["Go",go=True;While[go,psi2=Lh2[h1 . psi2];Pause[pause]; If[Abs[First[psi2]]>10*^-3,Break[]];n+=1],Method->"Queued"],Button["Stop",go=False,Method->"Preemptive"],Button["Reset",psi2=psi0;n=0;go=False;],Button["Update",{h1,Lh2,xx,psi0,vv}=prepareSolver[V,ψ0,xsolveDomain,δx,δt];psi2=psi0;go=False;n=0;]}}]},{Dynamic[If[ListQ[psi2],ListPlot[{joinCoordinates[xx,vv],joinCoordinates[xx,Re[psi2]],joinCoordinates[xx,Im[psi2]],joinCoordinates[xx,Abs[psi2]]},Joined-> True,PlotStyle-> {{Thick,Gray},{Thickness[0.001],Cyan},{Thickness[0.001],Pink},{Thickness[0.002],White}},Background-> Black,Axes-> False,PlotRange-> {{xvisDomain[[1]],xvisDomain[[2]]},{-1,2}},AspectRatio-> 3/4,ImageSize-> 800],Graphics[]],UpdateInterval-> pause]}}]