(* 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 cubic anharmonic oscillator x^2/2 + g^(1/2) x^3 ...\n"];
energy[0]=nqx+1/2;
Print["E[0] = ",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)]=<i|x^3|(i-3,i-1,i+1,i+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}];