Sample code 10.5 implements the variational Monte Carlo mothod to solve the 1-D TISE problem with a LJ potential, as suggested in Section 10.4, Chapter 10, Giodarno second edition. The expected GS energy is -1.8899 found using the shooting method.

Two type of difference equation for the TISE are used: (1) that using central differencing, and (2) that using single-step differencing. 

For central differencing method: It is found the algorithm is sufficiently accurate in producing ground state energy, EEstar = -2.07002 after 50000 step, NN=50, for the LJ potential with epsilon = 10, sigma = 1. 
 

For single-step method, result found are of poor quality. Not only the GS energy found are larger (Estar = -1.37713, 50000 timestep, NN=50), the shape of the GS solution is also found to be less satisfactory, see below.

(*Sample code 10.5*)
(* First, choose / define your potential *)
(* ONLY ONE of these options can be set to True *)
(* The rest must be set to False*)
varPotential1 = True;(* Harmonic potential *)
varPotential2 = False;
varPotential3 = False;
varPotential4 =False;
(*********** end of choosing potential type *)

If [varPotential1,(*Harmonic potential*)
k=1;xini=-2;xlast=2;
V[i_]:=0.5*k*x[i]^2;
];(*end if *)

If [varPotential2,(* LJ potential *)
xini=0.7;xlast=5.0;epsilon=10;sigma=1;
V[i_]:=4epsilon*((sigma/x[i])^12-(sigma/x[i])^6);
];(*end if *)
(* End of LJ potential *)


(*These are fixed parameters, need not change them unnecessarily*)
DeltaEE=10^8;EE=10^8;count=0;EEtrial=10^8;
x[i_]:=xini+i*Deltax;
(*End of These are fixed parameters, need not change them unnecessarily*)


(**** parameters, dont change it unnecessarily ***)
countlimit=2*10^4;(*This is the number of MC move to iterate. It is set to a large integer number. Default = 5*10^4. *)
NN=50;(*number of intervals*)
Deltax=(xlast-xini)/NN;
ilast=IntegerPart[(xlast-xini)/Deltax];m=1;
(**** parameters ***)

(*initialise phi[i]*)
For[i=-2,i≤ilast+2,i++,
phi[i]=0;
];(*end for i*)
i1=IntegerPart[(xini-xini)/Deltax];i4=IntegerPart[(xlast-xini)/Deltax];
For[i=i1,i≤i4,i++,
phi[i]=1;
];(*end for i*)

NC=Sum[Deltax*phi[i]^2,{i,0,ilast}];

(*EE=Sum[-(1/2)(1/m)*(phi[i]/Deltax)(phi[i+1]-2phi[i]+phi[i-1])+Deltax*V[i]*phi[i]^2,{i,0,ilast}];*)
EE=Sum[-(1/8)(1/m)*(phi[i]/Deltax)(phi[i+2]-2phi[i]+phi[i-2])+Deltax*V[i]*phi[i]^2,{i,0,ilast}];

EEstar=EE/NC;
txphi0=Table[{x[i],phi[i]},{i,0,ilast}];

(*Print["initially, Sum[phi[i],{i,0,ilast}]=",ilast*Deltax*Sum[phi[i]^2,{i,0,ilast}]];*)
(*end of initialist phi[i]*)

(* Keeping old phi before entering loop *)
(*For[i=0,i≤ilast,i++,
phiold[i]=phi[i]
];(*end for i*)*)
Print["initially, EEstar = ",EEstar];

While[count<countlimit,
count=count+1;

For[i=0,i≤ilast,i++,
phiold[i]=phi[i]
];(*end for i*)

iran=Random[Integer,{0,ilast}];
phiran=phi[iran];
phi[iran]=phiran+(2Random[Real]-1)(**Abs[phiran]*);

NC=Sum[Deltax*phi[i]^2,{i,0,ilast}];

(*EEtrial=Sum[-(1/2)(1/m)*(phi[i]/Deltax)(phi[i+1]-2phi[i]+phi[i-1])+Deltax*V[i]*phi[i]^2,
{i,0,ilast}];(*tested OK *)*)

EEtrial=Sum[-(1/8)(1/m)*(phi[i]/Deltax)(phi[i+2]-2phi[i]+phi[i-2])+Deltax*V[i]*phi[i]^2,{i,0,ilast}];(*tested, OK*)

EEstartrial=EEtrial/NC;

If[EEstartrial < EEstar,
EEstar=EEstartrial;
Print["Accept, count = ",{EEstar,count}];
,(*Else *)
(*restore phi to phiold*)
For[i=0,i≤ilast,i++,
phi[i]=phiold[i]
];(*end for i*)
];(*end if *)
(*************************************)
(*Print["count, NC, EEstartrial, EEstar, EEstartrial < EEstar ",{count, NC, EEstartrial, EEstar, EEstartrial < EEstar}];*)
(*NC=ilast*Deltax*Sum[phi[i]^2,{i,0,ilast}];
For[i=0,i≤ilast,i++,
phi[i]=phi[i]/Sqrt[NC];
];(*end for i*)
*)
];(*end while*)

NC=Sum[Deltax*phi[i]^2,{i,0,ilast}];
For[i=0,i≤ilast,i++,
phi[i]=phi[i]/Sqrt[NC];
];(*end for i*)
NC=Sum[Deltax*phi[i]^2,{i,0,ilast}];

(*EE=Sum[-(1/2)(1/m)*(phi[i]/Deltax)(phi[i+1]-2phi[i]+phi[i-1])+Deltax*V[i]*phi[i]^2,{i,0,ilast}];*)
EE=Sum[-(1/8)(1/m)*(phi[i]/Deltax)(phi[i+2]-2phi[i]+phi[i-2])+Deltax*V[i]*phi[i]^2,{i,0,ilast}];

EEstar=EE/NC;

Print["count,EEstar = ",{count,EEstar}];
txphi=Table[{x[i],phi[i]},{i,0,ilast}];
txV=Table[{x[i],V[i]},{i,i1,ilast}];
Print["ilast*Deltax*Sum[phi[i]^2,{i,0,ilast}]]=",ilast*Deltax*Sum[phi[i]^2,{i,0,ilast}]];
ListPlot[{txphi0,txphi,txV},Mesh->True,Joined->True,PlotRange->{{0,xlast},{-1,1}}];
ListPlot[{txphi,txV},PlotRange->{{xini-1,xlast+1},{-1,1}},Mesh->Automatic, Joined->True,PlotStyle->{PointSize->0.01}]