function [] = plot_edges_of_the_allowed_region_for_y_and_py(V, E, limits)

% region_edge_y_py(V, EM) plots the edges of the region of phase space
%   where V(0,y) is less than E, with y on the abscissa and the momentum on
%   the ordinate.  In other words, it plots the boundaries of the
%   allowed-energy region for the Poincare section of the phase space with
%   x=0 and px=0 (or, equivalently, x=0 and py=0, but it's more common to
%   plot y against its conjugate momentum).
%
%   The ranges for y and p are taken from the global variable
%   imDefaultLimits.  No, the code is not at all very generic. =)

m = 1;

NX = 50;  NY = 50;

Xstep = (limits(2)-limits(1))/(NX-1);
Xstep = limits(1):Xstep:limits(2);
Ystep = (limits(4)-limits(3))/(NY-1);
Ystep = limits(3):Ystep:limits(4);

for i = 1:NX,
  for j = 1:NY,
    data(j,i) = Ystep(j).^2./(2*m) + feval(V, 0, Xstep(i));
  end
end
contour(Xstep, Ystep, data, [E E]);  % must be specified twice to work
xlabel('y'); ylabel('p_y');  title(sprintf('E=%f', E));