(* 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["E<SUB>0</SUB> = ",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)]=<i|x|(i-4,i-2,i+2,i+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<SUB>",n,"</SUB> = ",energy[n]//InputForm];
   psi[[n+1,kaxh]]=0;
kbxmin0=kbxmin;
kbxmax0=kbxmax,{n,1,morder}];