(* Harmonic oscillator with omx=1, omy=2 on the lower surface *)
Clear["@"];
(********** INPUT starts **********)
{xmin,xmax}={-16.,16.};
{ymin,ymax}={-8.,8.};
{npx,npy}={482,241};
{nmin,nmax,nstep}={10,100,10};(* only even n! *)
tabdata={ (* x0, y0, omx, omy *)
   {0,0,0.02,0.18}, (* Fig. 1 *)
   {0,0,10.,2.2}, (* Fig. 2 *)
   {0,0,0.45,0.01}, (* Fig. 3 *)
   {0,0,2.,18.}, (* Fig. 4 *)
   {3.,0,2.,0.1}, (* Fig. 5 *)
   {3.,0,2.,10.} (* Fig. 6 *)
      };
(********** INPUT ends **********)
dirtemp="d:\\temp\\bilha";
sd=SetDirectory[dirtemp];
If[sd=!=dirtemp,Print["Directory ",dirtemp," does not exist. Exiting..."];
  Abort[]];
dx=(xmax-xmin)/npx;
dy=(ymax-ymin)/npy;
xpts=Table[xmin+(nx-1/2)dx,{nx,npx}];
ypts=Table[ymin+(ny-1/2)dy,{ny,npy}];
mtabdata=Length[tabdata];
psi[n_,x_,om_]:=(om/Pi)^(1/4)/(2^(n/2)Sqrt[n!])E^(-1/2 om x^2)*
    HermiteH[n,x Sqrt[om]];
psix=Table[psi[n,xpts[[nx]],1.],{n,0,nmax,2},{nx,npx}];
psiy=Table[psi[n,ypts[[ny]],2.],{n,0,nmax/2},{ny,npy}];
Do[
{x0,y0,omx,omy}=tabdata[[ntabdata]];
Print[StyleForm[
      StringForm[
"Fig. No. ``: \!\(x\_0\) = ``, \!\(y\_0\) = ``, \!\(\[Omega]\_x\) = ``, \!\(\[Omega]\_y\) = ``",
 ntabdata, Sequence @@ (#1 & ) /@ {x0, y0, omx, omy}], "Section"]];
Clear[int];
int[n_,om0_,om1_,a_]:=int[n,om0,om1,a]=
   E^(-a^2/2 om0 om1/(om0+om1))*
   om0^(1/4) om1^(1/4)/(om0+om1)^(n+1/2)*
   Sqrt[2^(n+1) n!]*
   Sum[2^(-2 i)/i! /(n-2 i)! (om1^2-om0^2)^i *
   If[n/2-i==0,1,(a^2 om0^2 om1)^(n/2-i)],{i,0,Floor[n/2]}];
Do[n=nn;
t=Timing[
zdata=Table[
   xn=xpts[[nx]];
   yn=ypts[[ny]];
   Sum[nqx=j;nqy=(n-j)/2;
      psix[[nqx/2+1,nx]]psiy[[nqy+1,ny]]*
      int[nqx,omx,1.,x0] int[nqy,omy,2.,y0]
   ,{j,0,n,2}],{ny,npy},{nx,npx}];
      ]//First;
Print["n = ",n," Time = ",t];
zmin0=Min[zdata];
zmax0=Max[zdata];
zscale=Max[Abs[zmin0],Abs[zmax0]];
pictfile="h2"<>ToString[ntabdata]<>"n"<>ToString[n]<>".tif";
cf=(z=zmin0+#(zmax0-zmin0);
Hue[If[z>0,0,0.65],1/4,z1=Abs[z/zscale]*1.1;If[z1<=1,z1,1]])&;
plt=ListDensityPlot[zdata, Mesh -> False,DisplayFunction->Identity,
PlotRange->{zmin0,zmax0},MeshRange->{{xmin,xmax},{ymin,ymax}},
   PlotLabel->
          StyleForm["n = "<>ToString[n],"Section",
            FontColor->RGBColor[1,0,0]],
   Background->RGBColor[1,1,0.9],ColorFunction->cf,AspectRatio->0.5];
Display[pictfile,plt,"TIFF",ImageSize->72 {8,4}];
,{nn,nmin,nmax,nstep}];
,{ntabdata,mtabdata}];