(* Input:  nqx,nqy - Oscillator quantum numbers
           omx,omy - Oscillator frequencies
           morder - Order of perturbation theory
           ndigits - Working precision
Example of input line: nqx=9; nqy=1; omx=1; omy=11/10; morder=12; ndigits=64;
*)

ncoef=morder+1;
energy[0]=(nqx+1/2)*omx+(nqy+1/2)*omy;
Print["E<SUB>0</SUB> = ",energy[0]//InputForm," (Unperturbed harmonic-oscillator energy)"];
   "(* PT coefficients for Barbanis potential *)";
mdouble=morder*2;
mbasisx=nqx+1*morder;
mbasisy=nqy+2*morder;
kax=nqx;
kay=nqy;
kayh=Floor[kay/2];
kymin=0;
If[kayh*2!=kay,kymin=1];
enxy=Table[0,{m,1,morder}];
xn=Table[Null,{kx,0,mbasisx},{k,-1,1,2}];
yn=Table[Null,{ky,kymin,mbasisy,2},{k,-2,2,2}];
kxm=0;
If[Floor[mbasisx/2]*2!=mbasisx,kxm=1];
psi=Table[Null,{n,0,mdouble},{kx,kxm,mbasisx,2},{ky,kymin,mbasisy,2}];
(* Calculating xn[i,(-1,1)]=<i|x|(i-1,i+1> and yn[i,(-2,0,2)]=<i|y^2|(i-2,i,i+2)> *)
Do[i1=i+1;
   xn[[i1,1]]=1;xn[[i1,2]]=(i+1)/(2*omx),{i,0,mbasisx}];
Do[ih=Floor[i/2]+1;
   yn[[ih,1]]=1;yn[[ih,2]]=(2*i+1)/(2*omy);
   yn[[ih,3]]=(i+1)*(i+2)/(4*omy^2),{i,kymin,mbasisy,2}];
(*     Harmonic oscillator wave function (zero order) *)
      psi[[1,Floor[nqx/2]+1,Floor[nqy/2]+1]]=1;
(*     Successive calculation of the expansion coefficients *)
kbxmin0=kax;
kbxmax0=kax;
kbymin0=kay;
kbymax0=kay;
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 *)
   If[nh*2==n,enxy[[nh]]=0];
   mdx=Min[n,mdouble-n];
   kbxmin=Max[kxmin,kax-mdx];
   kbxmax=kax+mdx;
   mdy=2*mdx;
   kbymin=Max[kymin,kay-mdy];
   kbymax=kay+mdy;
   Do[kbxh=Floor[kbx/2]+1;Do[a=0;kbyh=Floor[kby/2]+1;
ncmin=Max[Abs[kbx-kax],Ceiling[Abs[kby-kay]/2]];
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,kbyh]]*enxy[[nbch]],{nc,ncmin,n9,2}]];
kcxmin=Max[kbxmin0,kbx-1];
kcxmax=Min[kbxmax0,kbx+1];
kcymin=Max[kbymin0,kby-2];
kcymax=Min[kbymax0,kby+2];
Do[kcxh=Floor[kcx/2]+1;Do[kcyh=Floor[kcy/2]+1;
         a=a-psi[[n,kcxh,kcyh]]*xn[[kcx+1,(kbx-kcx+1)/2+1]]*yn[[kcyh,(kby-kcy)/2+2]],
         {kcy,kcymin,kcymax,2}],
      {kcx,kcxmin,kcxmax,2}];
If[kbx!=kax||kby!=kay,
        a=a/((kbx-kax)*omx+(kby-kay)*omy)];
psi[[n+1,kbxh,kbyh]]=a,{kby,kbymin,kbymax,2}],
      {kbx,kbxmin,kbxmax,2}];
   If[nh*2==n,
   enxy[[nh]]=-psi[[n+1,Floor[kax/2]+1,Floor[kay/2]+1]];
   energy[nh]=enxy[[nh]];
   Print["E<SUB>",nh,"</SUB> = ",energy[nh]//InputForm];
   If[enxy[[nh]]==Indeterminate,
      Print[];Print["Fermi resonance was encountered or insufficient precision. Higher order coefficients are indeterminate!"];Exit[]];
   psi[[n+1,Floor[kax/2]+1,Floor[kay/2]+1]]=0];
kbxmin0=kbxmin;
kbxmax0=kbxmax;
kbymin0=kbymin;
kbymax0=kbymax,{n,1,mdouble}];