function fig228
clc
clf
disp('Page 164 Fig2.28: Predator-Prey System')
disp('dR/dt = 2*R-1.2*R*F')
disp('dF/dt = -F+0.9*R*F')
disp('To be shown: vector fields (unnormalized)')
pause;

% Fig. 2.28
% Predator-Prey system
a = 2.;
b = -1.2;
c = -1.;
d = 0.9;
rmin = 0;
rmax = 4;
fmin = 0;
fmax = 4;
dr = 0.2;
df = 0.2;
[r,f] = meshgrid(rmin:dr:rmax,fmin:df:fmax);
u = a*r+b*r.*f;
v = c*f+d*r.*f;
leng = sqrt(u.^2+v.^2);
%
clf
% unnormalized version
quiver(r,f,u,v,0.5)
axis equal tight
xlabel 'r'
ylabel 'f'
title 'Vector Field for Predator-Prey system'
%
disp('To be shown: vector fields (normalized)')
pause
% normalize version
quiver(r,f,u./leng,v./leng,0.4)
axis equal tight
xlabel 'r'
ylabel 'f'
title 'Vector Field for Predator-Prey system'
%
disp('To be shown: equilibrium points')
pause
% solve for the equilibrium points
[x0,y0] = solve('2.*x0-1.2*x0*y0=0','-y0+0.9*x0*y0=0', 'x0','y0')
[x0 y0]
hold on;
plot(double(x0),double(y0),'go')
%
disp('To be shown: solution curves with initial condition (R,F)(0)=(1, 0.5)')
% solve for the solution curve
fode = inline('[2*y(1)-1.2*y(1)*y(2); -y(2)+0.9*y(1)*y(2)]', 't', 'y');
y10 = 1;
y20 = 0.5;
t0 = 0;
tmax = 20;
[t,ysol] = ode45(fode,[t0 tmax], [y10 y20]);
hold on
% solution curve in phase space
plot(ysol(:,1),ysol(:,2),'r-',double(x0),double(y0),'go')
hold off
%
pause;
subplot(1,2,1);
plot(ysol(:,1),ysol(:,2),'r-',double(x0),double(y0),'go')
subplot(1,2,2);
plot(t,ysol(:,1),'b-',t,ysol(:,2),'r-')
legend('r','f')
xlabel 't'
