%
%  code28.m -- gumdrop phase plane (djm - 26 nov 01)
%

%  coordinates
dx = 0.005;  dy = 0.005;
x = -1:dx:1;  y = -0.5+dy:dy:2.5;

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

first_int = (xg.^2) + yg - log(2*yg+1)/2;

%  2D phase plane

figure(1);  clf;
subplot(1,2,1)
cont = 0.02:0.1:1.02;
contour(x,y,first_int,cont)
hold on;  axis square;
xlabel('\bf u(t)')
ylabel('\bf u^\prime(t)')
title(['\bf gumdrop phase plane (first integral)'])
axis equal
axis([-1 1 -0.5 3.5])