(* 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 sextic anharmonic oscillator x^2/2 + g x^6 ...\n"];
energy[0]=nqx+1/2;
Print["E[0] = ",energy[0]//InputForm," (Unperturbed harmonic-oscillator energy)"];
mbasisx=nqx+3*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,-6,6,2}];
psi=Table[Null,{n,0,morder},{kx,kxmin,mbasisx,2}];
(* Calculating xn[i,(-6,-4,-2,0,2,4,6)]=<i|x|(i-6,i-4,i-2,i+2,i+4,i+6)> *)
Do[ih=(i-kxmin)/2+1;
xn[[ih,1]]=1;
xn[[ih,2]]=-9 + 6*i;
xn[[ih,3]]=15*(1 - i + i^2);
xn[[ih,4]]=5*(3 + 8*i + 6*i^2 + 4*i^3);
xn[[ih,5]]=15*(6 + 15*i + 14*i^2 + 6*i^3 + i^4);
xn[[ih,6]]=3*(120 + 298*i + 275*i^2 + 120*i^3 + 25*i^4 + 2*i^5);
xn[[ih,7]]=720 + 1764*i + 1624*i^2 + 735*i^3 + 175*i^4 + 21*i^5 + i^6,
   {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-6*mdx];
   kbxmax=kax+6*mdx;
   Do[a=0;kbxh=(kbx-kxmin)/2+1;
ncmin=Ceiling[Abs[kbx-kax]/6];
If[ncmin<=n9,
a=a+Sum[psi[[nc+1,kbxh]]*enx[[n-nc]],{nc,ncmin,n9}]];
kcxmin=Max[kbxmin0,kbx-6];
kcxmax=Min[kbxmax0,kbx+6];
Do[kcxh=(kcx-kxmin)/2+1;
a=a-psi[[n,kcxh]]*xn[[kcxh,(kbx-kcx)/2+4]],
      {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]]/8^n;
   Print["E[",n,"] = ",energy[n]//InputForm];
   psi[[n+1,kaxh]]=0;
kbxmin0=kbxmin;
kbxmax0=kbxmax,{n,1,morder}];