%% Predator-Prey Problem from Class

%% System of differential equations
% Let y1 denote the number of rabbits (prey), 
% let y2 denote the number of owls (predator).
% The population growth is described by
%
% * y1' = a1*y1
% * y2' = a2*y2
%
% Here the growth rate a1 for rabbits is positive for y2=0, but decreases
% with increasing y2;
% the growth rate a2 for owls is negative for y1=0, but increases with 
% increasing y1:
%
% * a1 = 2 - .5*y2
% * a2 = -1 + .5*y1
%
% Hence we obtain the following system of ODEs
%
% * y1' = (2-.5*y2)*y1
% * y2' = (-1+.5*y1)*y2

%% Given initial conditions
% The initial rabbit population is 6, the initial owl population is 2:
%
% * y1(0) = 6
% * y2(0) = 2

%% Definition of the function f
% We have an ODE system y' = f(t,y) where y is a vector and the vector-valued function f is given by 
f = @(t,y) [(2-.5*y(2))*y(1); (-1+.5*y(1))*y(2)];

%% Solve the initial value problem using ode45
% We want to find the functions y1(t) and y2(t) for t in [0,15]:
opt = odeset('RelTol',1e-10,'AbsTol',1e-10);  % use smaller tolerances
[ts,ys] = ode45(f,[0,15],[6;2],opt);

%% Plot both y1(t) and y2(t) vs. t
% It appears that the functions y1(t), y2(t) are periodic with a period
% of about 5.
plot(ts,ys)
legend('rabbits','owls'); xlabel('time t')

%% Plot y2 vs. y1 (Phase plane plot)
% The phase plane plot clearly shows that the solution is periodic
% since the trajectory is a closed curve.
%
% The vector field given by f shows the velocity vectors with which 
% the point (y1(t),y2(t)) moves along the trajectory. We see that
% the point moves along the closed curve _counterclockwise_ as t increases.
%
% (You first have to download the file
% <http://www.math.umd.edu/~tvp/246/vectfield.m vectfield.m>)
plot(ys(:,1),ys(:,2)); hold on
xlabel('rabbits'); ylabel('owls')
vectfield(f,0:.8:7.5,0:.8:10.5); hold off  % plot vector field


