(* Input:  nqx - Oscillator quantum number
           morder - Order of perturbation theory
Example of input line: nqx=0;morder=10;
*)
(* Testing that nqx and morder lie in a valid range *)
If[IntegerQ[nqx]&&nqx>=0,Print["Oscillator quantum number: n = ",nqx],
Print["Oscillator quantum number must be a non-negative integer!"];
Exit[]];
If[IntegerQ[morder]&&morder>0,Print["Order of perturbation theory: N = ",morder,"\n"],
Print["Order of perturbation theory must be a positive integer!"];
Exit[]];
ncoef=morder+1;
Print["Printing PT coefficients for the quartic anharmonic oscillator x^2/2 + g x^4 ...\n"];
energy[0]=nqx+1/2;
Print["E[0] = ",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[",n,"] = ",energy[n]//InputForm];
   psi[[n+1,kaxh]]=0;
kbxmin0=kbxmin;
kbxmax0=kbxmax,{n,1,morder}];