(* 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["E0 = ",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)]= and yn[i,(-2,0,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",nh," = ",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}];