(* Input: nqx - Oscillator quantum number morder - Order of perturbation theory Example of input line: nqx=0; morder=10; *) ncoef=morder+1; energy[0]=nqx+1/2; Print["E0 = ",energy[0]//InputForm," (Unperturbed harmonic-oscillator energy)"]; mbasisx=nqx+2*morder; kax=nqx; kaxh=Floor[kax/2]; kxmin=0; If[kaxh*2!=kax,kxmin=1]; enx=Table[0,{m,1,morder}]; xn=Table[Null,{kx,kxmin,mbasisx,2},{k,-4,4,2}]; psi=Table[Null,{n,0,morder},{kx,kxmin,mbasisx,2}]; (* Calculating xn[i,(-4,-2,0,2,4)]= *) Do[ih=(i-kxmin)/2+1; xn[[ih,1]]=1;xn[[ih,2]]=4*i-2; xn[[ih,3]]=3*(2*i^2+2*i+1); xn[[ih,4]]=(i+1)*(i+2)*(4*i+6); xn[[ih,5]]=(i+1)*(i+2)*(i+3)*(i+4),{i,kxmin,mbasisx,2}]; (* Harmonic oscillator wave function (zero order) *) kaxh=(kax-kxmin)/2+1; psi[[1,kaxh]]=1; (* Successive calculation of the expansion coefficients *) kbxmin0=kax; kbxmax0=kax; Do[n9=n-1; (* Calculating of wavefunctions at n-th step *) mdx=Min[n,morder-n]; kbxmin=Max[kxmin,kax-4*mdx]; kbxmax=kax+4*mdx; Do[a=0;kbxh=(kbx-kxmin)/2+1; ncmin=Ceiling[Abs[kbx-kax]/4]; If[ncmin<=n9, a=a+Sum[psi[[nc+1,kbxh]]*enx[[n-nc]],{nc,ncmin,n9}]]; kcxmin=Max[kbxmin0,kbx-4]; kcxmax=Min[kbxmax0,kbx+4]; Do[kcxh=(kcx-kxmin)/2+1; a=a-psi[[n,kcxh]]*xn[[kcxh,(kbx-kcx)/2+3]], {kcx,kcxmin,kcxmax,2}]; If[kbx!=kax, a=a/(kbx-kax)]; psi[[n+1,kbxh]]=a, {kbx,kbxmin,kbxmax,2}]; enx[[n]]=-psi[[n+1,kaxh]]; energy[n]=enx[[n]]/4^n; Print["E",n," = ",energy[n]//InputForm]; psi[[n+1,kaxh]]=0; kbxmin0=kbxmin; kbxmax0=kbxmax,{n,1,morder}];