(* ::Package:: *)


count=0;
Deltax=xlast/NN;
iini=IntegerPart[xini/Deltax];
ilast=IntegerPart[xlast/Deltax];x[0]=iota;
signphiprime=0;
DeltaE=DeltaEorig;




While[Abs[DeltaE]>tol,
For[i=iini,i<=ilast,i++,
xdum=i*Deltax;
V[i]=Potential[xdum]
];(*end for i*)

Enold=En;

(*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[2*phiL[i]-phiL[i-1]-2(En-V[i])(Deltax)^2*phiL[i]];
(*Print["i+1, phiL[i+1] ", i+1," ",phiL[i+1], " V[i]=",V[i]];*)

];(*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[2*phiR[i]-phiR[i+1]-2(En-V[i])(Deltax)^2*phiR[i]];
(*Print["i-1, phiR[i-1] ", i-1," xxx ",phiR[i-1]];*)
];(*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 both phiR and phiL  *)
(*scalefactorL=1/phiL[iMin];
scalefactorR=1/phiR[iMin];
For[i=ilast,i>=iini,i+=-1,
phiL[i]=phiL[i]*scalefactorL;
phiR[i]=phiR[i]*scalefactorR;
];(*end for i*)
(*Print[phiR[iMin]];
Print[phiL[iMin]];*)
(* End of doing the scalling on both phiR and phiL  *)
*)


(* Do the scalling on phiR only *)
scalefactor=phiL[iMin]/phiR[iMin];
For[i=ilast,i>=iini,i+=-1,
phiR[i]=phiR[i]*scalefactor
(*Print["phiR[i]=",phiR[i]];*)
];(*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];

If[Abs[DeltaE]<2*tol,
Print[
"En= ",En, " iMin=",iMin,
" iini= ",iini, " ilast = ",ilast,
" xini = ",xini, " xlast =",xlast,
"  DeltaE = ", DeltaE, " signphiprime=", signphiprime, 
" phiL[ilast]=", phiL[ilast], " phiR[iini]= ", phiR[iini]
];(*end print*)
];(*end if*)

(*Check that before accepting current solutino to DetalE, if 
phiL,R diverge at the edges, shift xini farther away from x=0 and redo the iteration 
by resetting DeltaE back to DeltaEorig*)
If[
Or[Abs[phiL[ilast]]> b , Abs[phiR[iini]]>b],
count=count+1;
DeltaE=DeltaEorig;
xini=xiniorig+0.1*count; iini=IntegerPart[xini/Deltax];
En=Enold;
xlast=xlastorig-0.1*count;ilast=IntegerPart[xlast/Deltax];
(*DeltaE=-DeltaE/2;*)
];(*end if*)


];(*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["n= ",n," 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}}]
*)
