%
%  duff_pp.m -- 17 jan 07 -- djm
%
%  - pendulum phase plane
%  - to run, type "code01" at the matlab prompt

global epsi
epsi = 0.2;

%  coordinates
N = 100;  L = 4;
the = -L:L/N:L;
phi = -L:L/N:L;

[Gthe Gphi] = meshgrid(the,phi);

%  first integral
Fint = (1/2)*Gphi.^2 + (1/2)*Gthe.^2 + (epsi/4)*Gthe.^4;

%  contourplot
figure(1);  clf;  hold on
contour(the,phi,Fint,[-2:1/2:6],'k')
axis equal;  axis([-1 1 -1 1]*L);

title('\bf duffing phase plane','fontsize',12)
xlabel('\bf y(t)-axis','fontsize',12)
ylabel('\bf y''(t)-axis','fontsize',12)