Page 130
% make an m-file named predPrey44.m with:
%   function Yprime=predPrey44(t,x)
%   r=1; a=1; m=1; b=1;
%   Yprime=[r*x(1)-a*x(1).*x(2);
%     -m*x(2)+b*x(1).*x(2)];
[t,Y]=ode23('predPrey44',[0 10],[1.5;0.5]); % ; for col vec
plot(t,Y) % both curves as the columns of Y vs t

               Page 131
x=Y(:,1); y=Y(:,2); plot(x,y)

               Page 134
% We must compute derivates numerically
% make an m-file predPrey44.m with:
%   function Yprime=predPrey44(t,x);
%   r=1; a=1; m=1; b=1;
%   Yprime=[r*x(1)-a*x(1).*x(2);
%     -m*x(2) + b*x(1).*x(2)];

% eps is matlab's smallest value; by divided diff.
% find the derivatives numerically; first at (0,0)
M1=(predPrey44(0,[eps 0]) - predPrey44(0,[0 0]))/eps;
   % this is the first column of the Jacobian at x=y=0
   % i.e. derivatives with respect to x
M2=(predPrey44(0,[0 eps]) - predPrey44(0,[0 0]))/eps;
   % the derivatives with respect to y
M=[M1 M2];  % the Jacobian
% calculate its eigenvalues
eig(M) % get 1 and -1, +1 means unstable at (0,0)

% now linearize at x=m/b=1, y=r/a=1
M1=(predPrey44(0,[1+eps 1])-predPrey44(0,[1 1]))/eps;
M2=(predPrey44(0,[1 1+eps])-predPrey44(0,[1 1]))/eps;
M=[M1 M2];
eig(M) % get +/-I (I=sqrt(-1)), so neutrally stable


                Exercises/Experiments

2.
% make an m-file, exer442.m
%   function Yprime=exer442(t,Y);
%   % Y(1)=x, Y(2)=y, Y(3)=z
%   Yprime=[-Y(1)+Y(1).*Y(2);
%     -Y(2)+2*Y(2).*Y(3)-Y(1).*Y(2);
%     2*Y(3)-Y(3).*Y(3)-Y(2).*Y(3)];

              For Grass
% grass
[t,Y]=ode23('exer442',[0 200],[0; 0; 1.5]);
plot(t,Y(:,3))

             For Grass and Sheep
% grass and sheep
[t,Y]=ode23('exer442',[0 200],[0; .5; 1.5]);
plot(t,Y)

             For Grass, Sheep, and Wolves
% all three
[t,Y]=ode23('exer442',[0 200], [.2; .5; 1.5]);
plot(t,Y(:,3),'g')  % grass behavior
hold on
plot(t,Y(:,2),'b')  % sheep behavior
plot(t,Y(:,1),'r')  % wolf behavior

5b.
% contents of the m-file exer445a.m:
%   function SIRprime=exer445a(t,SIR);
%   % S=SIR(1), I=SIR(2), R=SIR(3)
%   r=1; a=1;
%   SIRprime=[-r*SIR(1).*SIR(2);
%     r*SIR(1).*SIR(2)-a*SIR(2);a*SIR(2)];
r=1; a=1;
[t,SIR]=ode45('exer445a',[0 20], [2.8; .2; 0]);
plot(t,SIR)

5c.
N=100; % number steps per unit interval
delT=1/N; % so delta t=0.01
% t is now linked to index i by
% t=-1+(i-1)*delT where i=1,2,...,nFinal
% and the final index nFinal is given by
% solving tFinal = -1+(nFinal-1)*delT.
tFinal=5; nFinal=(tFinal+1)*N + 1;
% set up the initial values of R on -1 to 0
for i=1:N
   R(i)=0;
   S(i)=0; I(i)=0;
end
% work from t=0 in steps of delT
S(N+1)=2.8; I(N+1)=0.2; R(N+1)=0; 
for i=N+1:nFinal-1
   delY=delT*[-r*S(i)*I(i)+R(i-N);
      r*S(i)*I(i)-a*I(i); a*I(i)-R(i-N)];
   S(i+1) = S(i)+delY(1);
   I(i+1) = I(i)+delY(2);
   R(i+1) = R(i)+delY(3);
end 
% graph it
t=-1:delT:tFinal;
plot(t,S,t,I,t,R) % S blue, I green, R red