(* Calculation of pert. theory coefficients for arbitrary V(x) with a minimum *)
(* Input:  nqx - Oscillator quantum number
           morder - Order of perturbation theory;
           c[0], c[1], ..., c[2 morder+4] - coefficients of Taylor expansion of potential V(x) around its minimum
Example of input line: nqx=0; morder=10; c[0]=...;
*)

ncoef=morder;
nq=nqx;
ndig=128;
(
energy[0]=c[0];
Print["E<SUB>0</SUB> = ",energy[0]//pri," (Minimum of the potential)"];
If[morder==0,Goto[endanharm1]];
omega=Sqrt[2 c[2]];
energy[1]=(nqx+1/2)omega;
Print["E<SUB>1</SUB> = ",energy[1]//pri," (Quantum harmonic-oscillator energy)"];
If[morder==1,Goto[endanharm1]];
nm=ncoef-1;

(* R. - S. perturbation theory *)

na=nq+1;
one=1;
nmd=2*nm;
e=Table[0,{i,1,nm}];
f=Table[c[i+2]/omega^((i+2)/2),{i,1,nmd}];
nf=nmd;
two=2*one;
sq2=N[1/Sqrt[two],ndig];
np=Min[3*nm+1,na];
mm=np+3*nm;
mm2=mm+2;
mmh=Floor[(mm+1)/2];
npt=Floor[(np+1)/2];
t=Table[0,{n,1,nmd},{i,1,mmh}];
Do[
   n9=n-1;
   If[Abs[f[[n]]]>10.^-64,nf=n];
   in=1;
   If[Floor[(np+n)/2]*2==np+n,in=2];
   ma1=Min[np+3*n,np+3*(nmd-n)];
   mi1=Max[np-3*n,np-3*(nmd-n),in];
   Do[
      a=0;
      x=Table[0,{i,1,mm2+1}];
      nqt=Floor[(nq+1)/2];
      nb=na+nq-np;
      If[nq>1,x[[nq-1]]=(nb-one)*(nb-two)];
      x[[nq+1]]=2*nb-one;
      x[[nq+3]]=one;
      If[n==1,Goto[12]];
      Do[
         nt=n-nx;
         If[nx>nf,Goto[11]];
         in=1;
         If[Floor[(np+nt)/2]*2==np+nt,in=2];
         ma2=Min[np+3*nt,nq+nx+2];
         mi2=Max[np-3*nt,nq-nx-2,in];
         If[mi2>ma2,Goto[6]];
         c=0;
         Do[
            nrt=Floor[(nr+1)/2];
            b=x[[nr]]+(nr-np+na)*x[[nr+2]];
            x[[nr+1]]=b;
            c=c+b*t[[nt,nrt]], {nr,mi2,ma2,2}];
         a=a-c*f[[nx]]/2*sq2^nx;
         Label[11];nxh=Floor[nx/2];
         If[nxh*2!=nx,Goto[6]];
         a=a+e[[nxh]]*t[[nt,nqt]];
         Label[6];Null, {nx,1,n9}];
      Label[12];a=a-f[[n]]/2*sq2^n*(x[[np]]+na*x[[np+2]]);
      If[nq!=np,a=a/(nq-np)];
      t[[n,nqt]]=a, {nq,mi1,ma1,2}];
   nh=Floor[n/2];
   If[nh*2!=n,Goto[4]];
   e[[nh]]=-t[[n,npt]];
   energy[nh+1]=omega e[[nh]];
   Print["E<SUB>",nh+1,"</SUB> = ",energy[nh+1]//pri];
   t[[n,npt]]=0;
   Label[4];Null, {n,1,nmd}];
x=.;c=.;
Label[endanharm1];
);