%
%  code13.m -- cubic duffing phase plane (djm - 17 oct 01)
%

epsilon = [-0.25, 0.25];

%  parameters
Lx = 4.0;  Ly = 4.0;  Nx = 50;  Ny = 50;

%  coordinates
x = -Lx:Lx/Nx:Lx;  y = -Ly:Ly/Ny:Ly;

%  2D coordinates 
[xg,yg] = meshgrid(x,y);

for j=1:2
  first_int = (yg.^2 + xg.^2)/2 - epsilon(j)*(xg.^4)/4;
  cn = [-4:0.5:4];
  cont = (cn.^2)/2 - epsilon(j)*(cn.^4)/4;

%  2D plots

  figure(j);  clf;
  pcolor(x,y,first_int);  shading flat;  cbar
  hold on;  axis square
  contour(x,y,first_int,cont,'k')
  xlabel('\bf y')
  ylabel('\bf y^\prime')
  title(['\bf Duffing phase plane:  \epsilon=' num2str(epsilon(j))])
end