{x1, y1} = {x, y} - {x0, y0}; 
h = (x^2 + y^2)/2 + (h30*x^3 + h21*x^2*y + h12*x*y^2 + h03*y^3) + 
   (h40*x^4 + h31*x^3*y + h22*x^2*y^2 + h13*x*y^3 + h04*y^4); 
w = (omx^2*x1^2 + omy^2*y1^2)/2 + (w30*x1^3 + w21*x1^2*y1 + w12*x1*y1^2 + w03*y1^3) + 
   (w40*x1^4 + w31*x1^3*y1 + w22*x1^2*y1^2 + w13*x1*y1^3 + w04*y1^4); 
en = 1; 
{xmin, xmax, ymin, ymax} = {-2.1, 2.1, -2.1, 2.1}; 
dh = {D[h, x], D[h, y]}; 
dw = {D[w, x], D[w, y]}; 
s = N[NSolve[{Det[{dh, dw}] == 0, h == en}, {x, y}, WorkingPrecision -> 64]]; 
xy0 = xy = {}; 
ms = Length[s]; 
di = 10.^(-11); 
Do[sn = s[[n]];
   {xn, yn, wn} = {x, y, w} /. sn;
   {xnr, xni} = {Re[xn], Im[xn]}; 
   If[Abs[xni] > di, Continue[]];
   {ynr, yni} = {Re[yn], Im[yn]};
   If[Abs[yni] > di, Continue[]]; 
   xy = Append[xy, {xn, yn, wn}];
   If[xnr < xmin || xnr > xmax, Continue[]]; 
   If[ynr < ymin || ynr > ymax, Continue[]];
   xy0 = Append[xy0, {xn, yn, wn}];
  , {n, ms}]; 
xy = Sort[xy, #1[[3]] < #2[[3]] & ]; 
xy0 = Sort[xy0, #1[[3]] < #2[[3]] & ]; 
mxy = Length[xy]; mxy0 = Length[xy0];

Print["</P><P><PRE>Finding the constraint minimum ...</P>"];
If[mxy0==0,
   Print["<P><FONT COLOR=\"#FF0000\">No minimum found in the plotting area ",printF[xmin,5,2]," < x < ",printF[xmax,5,2],",     ",printF[ymin,5,2]," < y < ",printF[ymax,5,2],"</FONT>"];
   Print["<BR>List of all stationary points:"];
   Do[
   xymin=xy[[n]];
   Print["<P><I>x</I><SUB>",n,"</SUB> = ",printF[xymin[[1]],13,5],
     "       <I>y</I><SUB>",n,"</SUB> = ",printF[xymin[[2]],13,5],
     "       <I>W</I><SUB>",n,"</SUB> = ",printF[xymin[[3]],13,5] ]
      ,{n,mxy}];
   Print["<BR>Density plot aborted because (<I>x</I><SUB>min</SUB>,<I>y</I><SUB>min</SUB>) is out of plot range. Computations terminated.</PRE><HR>"];
   Exit[],
   xymin=xy0//First;
   Print["<P><I>x</I><SUB>min</SUB> = ",printF[xymin[[1]],8,5],
     "       <I>y</I><SUB>min</SUB> = ",printF[xymin[[2]],8,5],
     "       <I>W</I><SUB>min</SUB> = ",printF[xymin[[3]],8,5],
       "</PRE>"]
  ];