(* 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)"]; mdouble=morder*2; mbasisx=nqx+3*morder; kax=nqx; kaxmh=Floor[(kax+morder)/2]; kxmin=0; If[kaxmh*2!=kax+morder,kxmin=1]; enx=Table[0,{m,1,morder}]; xn=Table[Null,{kx,0,mbasisx},{k,-3,3,2}]; psi=Table[Null,{n,0,mdouble},{kx,kxmin,mbasisx,2}]; (* Calculating xn[i,(-3,-1,0,1,3)]= *) Do[i1=i+1;xn[[i1,1]]=1;xn[[i1,2]]=3*i; xn[[i1,3]]=3*(i+1)^2; xn[[i1,4]]=(i+1)*(i+2)*(i+3),{i,0,mbasisx}]; (* Harmonic oscillator wave function (zero order) *) kaxh=Floor[kax/2]+1; psi[[1,kaxh]]=1; (* Successive calculation of the expansion coefficients *) kbxmin0=kax; kbxmax0=kax; Do[n9=n-1; nh=Floor[n/2]; kanx=kax+n; kanxh=Floor[kanx/2]; kxmin=0; If[kanxh*2!=kanx,kxmin=1]; (* Calculating of wavefunctions at n-th step *) mdx=Min[n,mdouble-n]; kbxmin=Max[kxmin,kax-3*mdx]; kbxmax=kax+3*mdx; Do[a=0;kbxh=Floor[kbx/2]+1; ncmin=Ceiling[Abs[kbx-kax]/3]; nbcmax=n-ncmin; nbcmaxh=Floor[nbcmax/2]; If[nbcmaxh*2!=nbcmax,ncmin=ncmin+1]; If[ncmin<=n9, Do[nbc=n-nc; nbch=nbc/2; a=a+psi[[nc+1,kbxh]]*enx[[nbch]],{nc,ncmin,n9,2}]]; kcxmin=Max[kbxmin0,kbx-3]; kcxmax=Min[kbxmax0,kbx+3]; Do[kcxh=Floor[kcx/2]+1; a=a-psi[[n,kcxh]]*xn[[kcx+1,(kbx-kcx+1)/2+2]], {kcx,kcxmin,kcxmax,2}]; If[kbx!=kax, a=a/(kbx-kax)]; psi[[n+1,kbxh]]=a, {kbx,kbxmin,kbxmax,2}]; If[nh*2==n, enx[[nh]]=-psi[[n+1,kaxh]]; energy[nh]=enx[[nh]]/8^nh; Print["E",nh," = ",energy[nh]//InputForm]; psi[[n+1,kaxh]]=0]; kbxmin0=kbxmin; kbxmax0=kbxmax,{n,1,mdouble}];