(* 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 octic anharmonic oscillator x^2/2 + g x^8 ...\n"];
energy[0]=nqx+1/2;
Print["E[0] = ",energy[0]//InputForm," (Unperturbed harmonic-oscillator energy)"];
mbasisx=nqx+4*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,-8,8,2}];
psi=Table[Null,{n,0,morder},{kx,kxmin,mbasisx,2}];
(* Calculating xn[i,(-8,-6,-4,-2,0,2,4,6,8)]=<i|x|(i-8,i-6,i-4,i-2,i+2,i+4,i+6,i+8)> *)
Do[ih=(i-kxmin)/2+1;
xn[[ih,1]]=1;
xn[[ih,2]]=-20 + 8*i;
xn[[ih,3]]=14*(7 - 6*i + 2*i^2);
xn[[ih,4]]=28*(-3 + 7*i - 3*i^2 + 2*i^3);
xn[[ih,5]]=35*(3 + 8*i + 10*i^2 + 4*i^3 + 2*i^4);
xn[[ih,6]]=28*(30 + 83*i + 90*i^2 + 50*i^3 + 15*i^4 + 2*i^5);
xn[[ih,7]]=14*(360 + 990*i + 1073*i^2 + 600*i^3 + 185*i^4 + 30*i^5 +
     2*i^6);
xn[[ih,8]]=4*(5040 + 13788*i + 14896*i^2 + 8393*i^3 + 2695*i^4 +
     497*i^5 + 49*i^6 + 2*i^7);
xn[[ih,9]]=(1 + i)*(2 + i)*(3 + i)*(4 + i)*(5 + i)*(6 + i)*(7 + i)*
   (8 + i),
   {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-8*mdx];
   kbxmax=kax+8*mdx;
   Do[a=0;kbxh=(kbx-kxmin)/2+1;
ncmin=Ceiling[Abs[kbx-kax]/8];
If[ncmin<=n9,
a=a+Sum[psi[[nc+1,kbxh]]*enx[[n-nc]],{nc,ncmin,n9}]];
kcxmin=Max[kbxmin0,kbx-8];
kcxmax=Min[kbxmax0,kbx+8];
Do[kcxh=(kcx-kxmin)/2+1;
a=a-psi[[n,kcxh]]*xn[[kcxh,(kbx-kcx)/2+5]],
      {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]]/16^n;
   Print["E[",n,"] = ",energy[n]//InputForm];
   psi[[n+1,kaxh]]=0;
kbxmin0=kbxmin;
kbxmax0=kbxmax,{n,1,morder}];