Sample code 10.1. Shooting method to evaluate eigen values and eigen function of an 1D square well. The heigh of the well is V0.

tol=10^(-7);b=150;(*default = 150*)
DeltaE=0.1;V0=1000;(*default = 1000*)
NN=10000;a=1;Deltaa=0.5;(*default = 1*)
xini=-(1+Deltaa)a;xlast=(1+Deltaa)a;
Deltax=(2a)/NN;
Na=a/Deltax;
Nlimit=Na+IntegerPart[(xlast-a)/Deltax];
x[0]=0;
En=1.0;
(*Choose parity, either one*)
Parity=-1;
(*Parity=-1;*)
(*Choose parity, either one*)

If[Parity==1,phi[0]=1;phi[-1]=1];(*end of if Parity+1*)
If[Parity==-1,phi[0]=0;phi[-1]=-Deltax];(*end of if Parity+1*)


For[i=-Nlimit,i≤Nlimit,i++,
V[i]=0;
];(*end for i*)

For[i=Na+1,i≤Nlimit,i++,
V[i]=V0;
V[-i]=V0;
];(*end for i*)


lastdiverge=0;
While[Abs[DeltaE]>tol,
For[i=0,i≤Nlimit-1,i++,
x[i+1]=(i+1)*Deltax;
phi[i+1]=2phi[i]-phi[i-1]-2(En-V[i])(Deltax)^2*phi[i];
x[-(i+1)]=-x[i+1];
If[Abs[phi[i+1]]>b,Break[]];
];(*end for i*)
If[lastdiverge*Sign[phi[i]]<0,DeltaE=-DeltaE/2];(*end if*)
lastdiverge=Sign[phi[i]];
En=En+DeltaE;
Print["i = ",i,
" DeltaE ",DeltaE," ","lastdiverge= ",lastdiverge," ","phi[i] ",phi[i]];
];(*end while*)

(*for even parity case*)
If[Parity==1,
For[i=0,i≤Nlimit,i++,
phi[-(i)]=phi[i]
];(*end for i*)
];(*end if*)
(*for even parity case*)


(*for odd parity case*)
If[Parity==-1,
For[i=0,i≤Nlimit,i++,
phi[-(i)]=-phi[i]
];(*end for i*)
];(*end if*)

(*for odd parity case*)



Print["convergence achieved"]; 
Print["En= ",En];
txphiXN=Table[{N[x[i]],phi[i]},{i,-Nlimit,Nlimit}];
ListPlot[txphiXN,PlotLabel->"phi vs x, Not normlaised"]

(* The code below normalises the solution *)
Sumphi=2*Sum[phi[idum]^2*Deltax,{idum,-Nlimit,Nlimit}];
Print["Sumphi=",Sumphi];
For[i=-Nlimit,i≤Nlimit,i++,
phi[i]=phi[i]/Sqrt[Sumphi];
];(*end for i*)
NormalisedConst=Sum[phi[i]^2*Deltax,{i,-Nlimit,Nlimit}];
Print["NormalisedConst =",NormalisedConst];
txphiN=Table[{N[x[i]],phi[i]},{i,-Nlimit,Nlimit}];
ListPlot[txphiN,PlotLabel->"phi vs x. Normalised phi",PlotStyle->{Blue,Black}];
ListPlot[{txphiN,txphiXN},PlotLabel->"phi vs x. Unnormalised wv is Red",PlotStyle->{PointSize[0.0015],{Blue,Red}},PlotRange->{-2,2}]

Sample code 10.2 is modified from sample code 10.1 to evaluate eigen function phiL for Lenard-Jones Potential from the left. There is no need to impliment the shooting method here as phiL (phiR) does not need to obey the boundary condition at the far right (left) boundary. Evaluate phiL and phiR at the fixed En.

tol=10^(-7);b=10;sigma=1;epsilon=10;
DeltaE=0.1;
NN=10000;
xini=0.7;(*default = 0.7*)
xlast=5;
Deltax=xlast/NN;
iini=IntegerPart[xini/Deltax];
ilast=IntegerPart[xlast/Deltax];

x[0]=iota;iota=10^(-8);
En=-1.6;
;
Potential[xdum_]:=4epsilon*((sigma/xdum)^12-(sigma/xdum)^6);

V[0]=Potential[iota];
For[i=1,i≤ilast,i++,
xdum=i*Deltax;
V[i]=Potential[xdum]
];(*end for i*)

(*Initial condition for phiL*)
phiL[iini]=0;phiL[iini-1]=phiL[iini]-Deltax*iota;
For[i=iini,i≤ilast-1,i++,
phiL[i+1]=2phiL[i]-phiL[i-1]-2(En-V[i])(Deltax)^2*phiL[i];
(*Print[i+1," ",phiL[i+1]];*)
];(*end for i*)


(*Initial condition for phiR*)(* Count from the Right *)
phiR[ilast]=0;phiR[ilast-1]=phiR[ilast]-Deltax*iota;
For[i=ilast-1,i≥iini,i+=-1,
phiR[i-1]=2phiR[i]-phiR[i+1]-2(En-V[i])(Deltax)^2*phiR[i];
];(*end for i*)

Print["En= ",En];

txphiL=Table[{N[xini+(i-1)*Deltax],phiL[i]},{i,iini,ilast}];
txphiR=Table[{N[xlast+(i-1)*Deltax],phiR[i]},{i,1,ilast}];

ListPlot[{txphiL,txphiR},PlotLabel->"phiL,phiR vs x, Not normlaised. phiL is Blue",PlotStyle->{PointSize[0.0015],{Blue,Red}},PlotRange->All]

ListPlot[{txphiL},PlotLabel->"phiL vs x, Not normlaised.",PlotStyle->{PointSize[0.0015]}];
ListPlot[{txphiR},PlotLabel->"phiR vs x, Not normlaised.",PlotStyle->{PointSize[0.0015]}];

Sample code 10.3: Full implimentation of the matching method for the TISE with a LJ potential.

tol=10^(-7);b=10;sigma=1;epsilon=10;
NN=7000;(*default 10000*)
xini=0.65;(*Default 0.65*)
xlast=5;b=100;
Deltax=xlast/NN;
iini=IntegerPart[xini/Deltax];
ilast=IntegerPart[xlast/Deltax];

x[0]=iota;iota=10^(-4);
En=-2;signphiprime=0;DeltaE=0.1;
;
Potential[xdum_]:=4epsilon*((sigma/xdum)^12-(sigma/xdum)^6);

V[0]=Potential[iota];
For[i=1,i≤ilast,i++,
xdum=i*Deltax;
V[i]=Potential[xdum]
];(*end for i*)


(* Find xmin, x where V is minimum *)
xmin=Minimize[Potential[x],x>0,x];
(*xMin=x/.xmin[[2]][[1]];*)
xMin=1.0;
iMin=IntegerPart[xMin/Deltax];
Plot[Potential[x],{x,xMin-0.1,xMin+0.1}];


While[Abs[DeltaE]>tol,
(*Initial condition for phiL*)
phiL[iini]=0;phiL[iini-1]=phiL[iini]-Deltax*iota;
For[i=iini,i≤ilast-1,i++,
phiL[i+1]=N[2phiL[i]-phiL[i-1]-2(En-V[i])(Deltax)^2*phiL[i]];
(*If[Abs[phiL[i+1]]>b,
Break[]];*)
(*Print[i+1," ",phiL[i+1]];*)
];(*end for i*)

(*Initial condition for phiR*)(* Count from the Right *)
phiR[ilast]=0;phiR[ilast-1]=phiR[ilast]-Deltax*iota;
For[i=ilast-1,i≥iini,i+=-1,
phiR[i-1]=N[2phiR[i]-phiR[i+1]-2(En-V[i])(Deltax)^2*phiR[i]];
(*If[Abs[phiR[i-1]]>b,Break[]];*)
];(*end for i*)

txphiL:=Table[{N[xini+(i-1)*Deltax],phiL[i]},{i,iini,ilast}];
txphiR:=Table[{N[xlast+(i-1)*Deltax],phiR[i]},{i,1,ilast}];

ListPlot[{txphiL,txphiR},PlotLabel->"phiL,phiR vs x, Not normlaised. phiL is Blue",PlotStyle->{PointSize[0.0015],{Blue,Red}}];

ListPlot[{txphiL},PlotLabel->"phiL vs x, Not normlaised.",PlotStyle->{PointSize[0.0015]}];
ListPlot[{txphiR},PlotLabel->"phiR vs x, Not normlaised.",PlotStyle->{PointSize[0.0015]}];
(*****************************************************************)



(* Do the scalling on phiR only *)
scalefactor=phiL[iMin]/phiR[iMin];
For[i=ilast,i≥iini,i+=-1,
phiR[i]=phiR[i]*scalefactor
];(*end for i*)
(*Print[phiR[iMin]];
Print[phiL[iMin]];*)
(* End of doing the scalling on phiR only *)


(* Now check the matching of phi'[ixMin]*)
phiLprime=(phiL[iMin+1]-phiL[iMin-1])/(2Deltax);
phiRprime=(phiR[iMin+1]-phiL[iMin-1])/(2Deltax);
(*Print["iMin=",iMin];
Print["phiLprime = ",phiLprime, " phiRprime=", phiRprime];*)
If[signphiprime*Sign[phiLprime-phiRprime]<0,DeltaE=-DeltaE/2];(*end if*)
En=En+DeltaE;
signphiprime=Sign[phiLprime-phiRprime];
Print["En= ",En, "  ","DeltaE = ",DeltaE, "   signphiprime=",signphiprime, " phiL[ilast]=",phiL[ilast], " phiR[iini]= ",phiR[iini]
];
];(*end while*)


(* The code below normalises the solution *)
SumphiL=Sum[phiL[i]^2*Deltax,{i,iini,ilast}];
Print["SumphiL=",SumphiL];
SumphiR=Sum[phiR[i]^2*Deltax,{i,iini,ilast}];
Print["SumphiR=",SumphiR];

For[i=iini,i≤ilast,i++,
phiL[i]=phiL[i]/Sqrt[SumphiL];
phiR[i]=phiR[i]/Sqrt[SumphiR];
];(*end for i*)
NormalisedConstL=Sum[phiL[i]^2*Deltax,{i,iini,ilast}];
Print["NormalisedConstL =",NormalisedConstL];
NormalisedConstR=Sum[phiR[i]^2*Deltax,{i,iini,ilast}];
Print["NormalisedConstR =",NormalisedConstR];
Print["En=",En];
Print["xMin=",N[xMin]];
txphiL=Table[{N[i*Deltax],phiL[i]},{i,iini,ilast}];
txphiR=Table[{N[i*Deltax],phiR[i]},{i,iini,ilast}];
txV=Table[{i*Deltax,N[V[i]]},{i,iini,ilast}];

ListPlot[{txphiL,txphiR,txV},PlotLabel->"phiL,phiR, V vs x, Normlaised. phiL is Green",PlotStyle->{PointSize[0.0015],Green,Red,Black},PlotRange->{{0,5},{-10,3}}]
ListPlot[{txphiL,txV},PlotLabel->"phiL vs x, Normlaised.",PlotStyle->{PointSize[0.0015]},PlotRange->{{0,5},{-10,3}}]
ListPlot[{txphiR,txV},PlotLabel->"phiR vs x, Normlaised.",PlotStyle->{PointSize[0.0015]},PlotRange->{{0,5},{-10,3}}]
ListPlot[{txphiL,txphiR},PlotLabel->"phiL,phiR, vs x, Normlaised. phiL is Green",PlotStyle->{PointSize[0.0015],Green,Red,Black},PlotRange->{{0,5},{-1,3}}]