***eigenfunctions***

psinodd[a_,n_,x_]:=Sqrt[2/a]Cos[n Pi x/a]

psineven[a_,n_,x_]:=Sqrt[2/a]Sin[n Pi x/a]

psi[a_,n_,x_]:=If[OddQ[n],psinodd[1,n,x],psineven[1,n,x]]

Manipulate[Plot[(psi[1,n,x])^2,{x,-1/2,1/2},PlotStyle->{Red,Thick},AxesLabel->{Style["x/a",Large],Style["a|ψ(x/a)((|)^2)",Large]},LabelStyle->{FontSize->16}],{n,1,50}]

Manipulate[Plot[2 Sin[n Pi x]^2 ,{x,0,1},PlotStyle->{Red,Thick},AxesLabel->{Style["x/a",Large],Style["a|(Subscript[ψ, n])(x/a)((|)^2)",Large]},LabelStyle->{FontSize->16},PlotRange->{0,2.2}],{n,1,20}]

***momentum eigenfunctions***

Integrate[Sqrt[1/(Pi 2 )]Sqrt[2/a]Sin[n Pi x/a] Exp[I k x],{x,0,a}]

phi[a_,n_,k_]:=((-n π+(E^(I a k)) (n π Cos[n π]-I a k Sin[n π]))/((Sqrt[(1/a)]) (Sqrt[π]) (a k-n π) (a k+n π)))

Manipulate[Plot[Abs[phi[1,n,k]]^2,{k,- n Pi-6,n Pi+6},PlotStyle->{Red,Thick},AxesLabel->{Style["ka",Large],Style["\[LeftBracketingBar](Subscript[ϕ, n])(ka)(\[RightBracketingBar]^2)/a",Large]},LabelStyle->{FontSize->16},PlotRange->All],{n,1,20}]

***Δx  and Δk***

xprob=(2/a) Sin[n Pi x/a]^2

xav=Integrate[xprob x,{x,0,a}]

Simplify[%,n∈Integers]

Δx=Sqrt[Integrate[xprob (x-a/2)^2,{x,0,a}]]

Simplify[%,n∈Integers]

uncx[n_]:=a((Sqrt[(( (-6+(n^2) (π^2)))/(n^2))])/(2 (Sqrt[3]) π))

Simplify[((-n π+(E^(I a k)) (n π Cos[n π]-I a k Sin[n π]))/((Sqrt[(1/a)]) (Sqrt[π]) (a k-n π) (a k+n π)))((-n π+(E^(-I a k)) (n π Cos[n π]+I a k Sin[n π]))/((Sqrt[(1/a)]) (Sqrt[π]) (a k-n π) (a k+n π))),Assumptions->n∈Integers]

2((a  (1-(((-1))^n) Cos[a k]) (n^2) π)/(((((a k)^2-(n π)^2))^2) ))

kprob[n_,a_,k_]:=((2 a (n^2) π (1-(((-1))^n) Cos[a k]))/((((a^2) (k^2)-(n^2) (π^2)))^2))

Limit[kprob[n,a,k],k->n Pi/a,Assumptions->n∈Integers]

Δk1=Integrate[kprob[1] k^2,{k,-Infinity,Infinity},Assumptions->a>0]

Δk2=Integrate[kprob[2] k^2,{k,-Infinity,Infinity},Assumptions->a>0]

***note that Δk=nℼ/a since it is related to square root of energy***

unck[n_]:=n Pi/a

unc[n_]:=uncx[n] unck[n]

Simplify[unc[n],n>0]

unc[n_]:=(1/2) (Sqrt[-2+(((n^2) (π^2))/3)])

unc[1.]

unc[2.]

***evolution of a wave packet***

***take psi0=nor[a,b]Exp[-(x-a/2)^2/(2 b^2)] Exp[I k0 x]****

Integrate[Exp[-(x-a/2)^2/( b^2)] ,{x,0,a}]

nor[a_,b_]:=1/Sqrt[b (Sqrt[π]) Erf[(a/(2 b))]]

psi0[a_,b_,k0_,x_]:=nor[a,b]Exp[-(x-a/2)^2/( 2b^2)] Exp[I k0 x]

**get expansion coefficients***

Integrate[Sqrt[(2/a)] Sin[n Pi x/a]psi0[a,b,k0,x],{x,0,a}]

c[a_,b_,n_,k0_]:=(1/(2 (Sqrt[b Erf[(a/(2 b))]])))(Sqrt[(1/a)]) b (E^(-((-I (a^3) k0+2 a (b^2) k0 n π+(b^2) (n^2) (π^2)+(a^2) ((b^2) (k0^2)+I n π))/(2 (a^2))))) (π^(1/4)) (-(E^((2 (b^2) k0 n π)/a)) Erfi[((-I (a^2)+2 a (b^2) k0-2 (b^2) n π)/(2 (Sqrt[2]) a b))]+(E^((2 (b^2) k0 n π)/a)) Erfi[((I (a^2)+2 a (b^2) k0-2 (b^2) n π)/(2 (Sqrt[2]) a b))]+(E^(I n π)) (Erfi[((-I (a^2)+2 a (b^2) k0+2 (b^2) n π)/(2 (Sqrt[2]) a b))]-Erfi[((I (a^2)+2 a (b^2) k0+2 (b^2) n π)/(2 (Sqrt[2]) a b))]))

psit[a_,b_,x_,t_,k0_,nmin_,nmax_]:=Sum[c[a,b,n,k0] Sqrt[(2/a)] Sin[n Pi x/a]Exp[-2 Pi I n^2 t],{n,nmin,nmax}]

***t is measured in units of 2 m a^2/(π^2)/2*Pi****

Plot[{Abs[psit[1,.1,x,0,0,1,3]]^2,Abs[psi0[1,.1,0,x]]^2},{x,0,1},PlotRange->All]

Plot[{Abs[psit[1,.1,x,0,0,1,15]]^2,Abs[psi0[1,.1,0,x]]^2},{x,0,1},PlotRange->All]

***with 15 terms, no difference, the narrower the packet, the more terms needed****

Plot[{Abs[psit[1,.1,x,0,20 Pi,1,15]]^2,Abs[psi0[1,.1,0,x]]^2},{x,0,1},PlotRange->All]

Plot[{Abs[psit[1,.1,x,0,20 Pi,5,35]]^2,Abs[psi0[1,.1,0,x]]^2},{x,0,1},PlotRange->All]

***If (Subscript[k, 0])=20π and a=1, must take n=20±15 for convergence***

Plot[{Abs[psit[1,.1,x,.02,5 Pi,0,20]]^2,Abs[psi0[1,.1,0,x]]^2},{x,0,1},PlotRange->All]

***approximate evolution of a wave packet***

***take psi0=nor[a,b]Exp[-(x-a/2)^2/(2 b^2)] Exp[I k0 x]****

Integrate[Exp[-(x-a/2)^2/( b^2)] ,{x,0,a}]

nor[a_,b_]:=1/Sqrt[b (Sqrt[π]) Erf[(a/(2 b))]]

psi0[a_,b_,k0_,x_]:=nor[a,b]Exp[-(x-a/2)^2/( 2b^2)] Exp[I k0 x]

**get expansion coefficients***

Integrate[Sqrt[(2/a)] Sin[n Pi x/a]psi0[a,b,k0,x],{x,-Infinity,Infinity},Assumptions->k0>0&&n>0&&a>0&&b>0&&a>b]

capp[a_,b_,n_,k0_]:=((I b (E^(-((-I (a^3) k0+2 a (b^2) k0 n π+(b^2) (n^2) (π^2)+(a^2) ((b^2) (k0^2)+I n π))/(2 (a^2))))) (-(E^(I n π))+(E^((2 (b^2) k0 n π)/a))) (π^(1/4)))/(Sqrt[a b Erf[(a/(2 b))]]))

capp[1,.1,3,.73]

psitapp[a_,b_,x_,t_,k0_,nmin_,nmax_]:=Sum[capp[a,b,n,k0] Sqrt[(2/a)] Sin[n Pi x/a]Exp[-2 Pi I n^2 t],{n,nmin,nmax}]

***t is measured in units of 2 m a^2/(π^2)****

Plot[{Abs[psitapp[1,.1,x,0,0,1,3]]^2,Abs[psi0[1,.1,0,x]]^2},{x,0,1},PlotRange->All]

Plot[{Abs[psitapp[1,.1,x,0,0,1,15]]^2,Abs[psi0[1,.1,0,x]]^2},{x,0,1},PlotRange->All]

***with 15 terms, no difference, the narrower the packet, the more terms needed****

Plot[{Abs[psitapp[1,.1,x,0,20 Pi,1,15]]^2,Abs[psi0[1,.1,0,x]]^2},{x,0,1},PlotRange->All]

Plot[{Abs[psitapp[1,.1,x,0,20 Pi,5,35]]^2,Abs[psi0[1,.1,0,x]]^2},{x,0,1},PlotRange->All]

***If (Subscript[k, 0])=20π and a=1, must take n=20±15 for convergence***

Manipulate[Plot[{Abs[psitapp[1,.1,x,t,5 Pi,0,20]]^2,Abs[psi0[1,.1,0,x]]^2},{x,0,1},PlotRange->All],{t,0,1}]

****at origin****

psi0sym[a_,b_,k0_,x_]:=(b (Sqrt[π]))^(-1/2)Exp[-(x)^2/( 2b^2)] Exp[I k0 x]

psi0sym[1,.1,0,0]

Integrate[Exp[-x^2/b^2] ,{x,-Infinity,Infinity},Assumptions->b>0]

**get expansion coefficients***

Integrate[Sqrt[(2/a)] Cos[(2 n-1) Pi x/a]psi0sym[a,b,k0,x],{x,-Infinity,Infinity},Assumptions->k0>0&&n>0&&a>0&&b>0&&a>b]

ansym[a_,b_,n_,k0_]:=((b (E^(-(b^2) ((k0^2)+(((((1-2 n))^2) (π^2))/(a^2))))) ((E^(((b^2) (((a k0+π-2 n π))^2))/(2 (a^2))))+(E^(((b^2) (((a k0+(-1+2 n) π))^2))/(2 (a^2))))) (π^(1/4)))/(Sqrt[a b]))

Integrate[Sqrt[(2/a)] Sin[(2 n) Pi x/a]psi0sym[a,b,k0,x],{x,-Infinity,Infinity},Assumptions->k0>0&&n>0&&a>0&&b>0&&a>b]

bnsym[a_,b_,n_,k0_]:=((I b (E^(-(((b^2) (((a k0+2 n π))^2))/(2 (a^2))))) (-1+(E^((4 (b^2) k0 n π)/a))) (π^(1/4)))/(Sqrt[a b]))

bnsym[1,.1,1,0]

psitsym[a_,b_,x_,t_,k0_,nmin_,nmax_]:=Sum[ansym[a,b,n,k0] Sqrt[(2/a)] Cos[(2 n-1) Pi x/a]Exp[-2 Pi I (2 n-1)^2 t]+bnsym[a,b,n,k0] Sqrt[(2/a)] Sin[2n Pi x/a]Exp[-2 Pi I (2 n)^2 t],{n,nmin,nmax}]

Manipulate[Plot[{Abs[psitsym[1,.1,x,t,0,1,15]]^2,Abs[psi0sym[1,.1,0,x]]^2},{x,-1/2,1/2},PlotRange->All],{t,0,1}]

Manipulate[Plot[{Abs[psitsym[1,.1,x,t,5 Pi,1,20]]^2,Abs[psi0sym[1,.1,0,x]]^2},{x,-1/2,1/2},PlotRange->All],{t,0,1}]