Functions

MakePsi0[tildep_, tildeb_,Nn_] := Module[{ret,dq},
dq = 1./(Nn+1);
ret = Table[(2 tildeb/Pi)^(1/4) Exp[-tildeb ( j dq - 1/2)^2] Exp[I tildep (j dq-1/2)],{j,1,Nn}];
Return[ret];
]

PlotPsi[inPsisteps_, inmovsteps_] := Module[{retplots},
retplots =Table[ListLinePlot[ Abs[inPsisteps[[j]]]^2,PlotRange->{0,10}],{j,1,Length[inPsisteps],inmovsteps}];
Return[retplots];
]

MakeH[Vtilde_, Nn_] := Module[{Hmat,dq,index},
dq = 1/(Nn + 1);
Hmat = Table[0.0,{i,1,Nn},{j,1,Nn}];
For[index = 2, index ≤ Nn-1, index=index+1,
Hmat[[index, index-1;;index+1]] = {-1/dq^2,(+2/dq^2 + Vtilde[index dq ] ), -1/dq^2};
];
Hmat[[1,1;;2]] = {2/dq^2 +  Vtilde[dq] ,-1/dq^2};
Hmat[[Nn,Nn-1;;Nn]] = {-1/dq^2,2/dq^2 + Vtilde[Nn dq] };
Return[Hmat];
]

RunCN[ds_, Nn_,Ns_, tildeb_,tildep_,tildeV_] := Module[{ Psiout,Psiin, Prop,index,H,Q,II,Pbar,dq},
dq= 1./(Nn+1);
Psiin = MakePsi0[tildep, tildeb,Nn];
H = MakeH[tildeV, Nn];
II = DiagonalMatrix[Table[1.0,{j,1,Length[H]}]];
Pbar = II + I ds H/2.;
Prop =(Inverse[Pbar].Conjugate[Pbar]);
Psiout = {Psiin};
For[index = 1, index ≤ Ns, index = index+1,
Psiin = Prop.Psiin;
Psiout = Join[Psiout,{Psiin}];
];
Return[Psiout];
]

CheckNorm[inPsisteps_] := Module[{retnorm,dq},
dq = 1./(Length[inPsisteps[[1]]]+1);
retnorm = Table[inPsisteps[[jj]].Conjugate[inPsisteps[[jj]]] dq,{jj,1,Length[inPsisteps]}];
Return[retnorm];
]

No Potential

Free Spreading

Psi0 = MakePsi0[0., 80,100];

ListPlot[Psi0]

V[x_]:= 0

test = RunCN[.0001,100,100,80, 0,V];

ListPlot[CheckNorm[test]]

mov = PlotPsi[test,1];

ListAnimate[mov]

Moving Right

Psi0 = MakePsi0[0., 80,100];

ListPlot[Psi0]

V[x_]:= 0

test = RunCN[.0001,100,100,80, 10,V];

ListPlot[CheckNorm[test]]

mov = PlotPsi[test,1];

ListAnimate[mov]

Gravity

Drop from Rest

V[x_]:= 10000x

test = RunCN[.0001,100,400,80, 0,V];

mov = PlotPsi[test,1];

ListAnimate[mov]

ListPlot[CheckNorm[test]]

Thrown upward

V[x_]:= 10000x

test = RunCN[.0001,100,200,80, 40,V];

mov = PlotPsi[test,1];

ListAnimate[mov]

ListPlot[CheckNorm[test]]