(* ::Package:: *)

tol=10^(-5); count=0;
NN=5000;(*default 10000*)
xiniorig=0.7;(*Default 0.86*)
DeltaEorig=0.1;
b=100;

xini=xiniorig;
xlast=5; 
Deltax=xlast/NN;
iini=IntegerPart[xini/Deltax];
ilast=IntegerPart[xlast/Deltax];

x[0]=iota;iota=10^(-4);
(*En=-2*)
signphiprime=0;

DeltaE=DeltaEorig;

<<potential.m;




(* 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,

For[i=iini,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]=N[2phiL[i]-phiL[i-1]-2(En-V[i])(Deltax)^2*phiL[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[2phiR[i]-phiR[i+1]-2(En-V[i])(Deltax)^2*phiR[i]];
];(*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
];(*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, " xini = ",xini, "  ","DeltaE = ",DeltaE, 
"   signphiprime=",signphiprime, " phiL[ilast]=",phiL[ilast], " phiR[iini]= ",phiR[iini]];
];(*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;
xini=xiniorig+0.05*count; iini=IntegerPart[xini/Deltax];DeltaE=DeltaEorig;
];(*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}}]
*)
