% This file contains MATLAB code to guide you through problem set 2
% and more generally to teach you phase plane diagram analysis for
% differential equation.  Try changing the parameters and see how 
% these changes effect the solutions.

clf;

global alpha beta;
more on; echo on;

% Consider the following system of first order differential equations
%
%  dy_1/dt = (1 - alpha y_2) y_1
%
%  dy_2/dt = (- 1 + beta y_1) y_2
%
%  where, variables y_1 and y_2 measure the sizes of the prey and 
%  predator populations, respectively.  The quadratic cross terms 
%  account for the interaction between the species.  This system is
%  known as the LOTKA-VOLTERRA predator-prey model.

% The first task when analyzing the system of differential equations
% is to plot both variables as a function of time.  For this one needs to
% specify an interval of time, initial conditions of the size of the prey
% and predator population and values for the parameters alpha and
% beta and then one needs to solve the system. 

echo off;

T = input('What time interval do you want for the simulation? T=');
y_1_init = input('Specify an initial size of the predator population y_1_init=');
y_2_init = input('Specify an initial size of the prey population y_2_init=');
alpha = input('Pick a value for the parameter alpha alpha=');
beta = input('Pick a value for the parameter beta beta=');

echo on;

% To solve the differential equations we use RUNGE KUTTA method referred 
% to as ode23. To read more about this code type help ode23 inside MATLAB.
% Function "prey" defines the differential equations, tspan is the period
% of simulation and y0 are the initial population values.  ode23 returns two
% arguments: t - different points in time between 0 and T; and 
% y - value of the two populations at those times;

y0 = [y_1_init; y_2_init];  % Initial Values;
tspan = [0 T];              % Simulation Period; 

[t y] = ode23('prey',tspan,y0);

% We can now plot the time paths of the two populations.

plot(t,y)

title('Lotke-Volterra Equation Time History');
xlabel('TIME');
ylabel('POPULATIONS');
legend(['y_1';'y_2']);

pause;

% The next task is to undertake the phase-plane analysis.  For this we first
% plot the isoclines.  The isoclines are obtained by setting the differential
% equations to 0.  Thus the isocline dy_1/dt = 0 is given by y_2 = 1/alpha
% while the isocline dy_2/dt = 0 is given by y_1 = 1/beta;

clf;

% Plot of the dy_1/dt = 0 isocline;

plot([0; round(max(y(:,1)))], [1/alpha;1/alpha]);

hold on;

% Plot of the dy_2/dt = 0 isocline;

plot([1/beta;1/beta],[0;round(max(y(:,2)))]);

xlabel('y_1');
ylabel('y_2');

% Next we determine the direction of the trajectory in each isosector
% (quadrant formed by the isoclines). For this we pick a point in each
% quadrant and determine the sign the differential equation at that point.  A
% negative value implies that the variable is decreasing while a positive
% value implies that the variable is increasing in that isosector.

% For the north-east isosector:

echo off;

p_1 = input('Pick a point on the horizontal axis in the isosector p_1=');
p_2 = input('Pick a point on the vertical axis in the isosector p_2=');

echo on;

   % The directions of the variables are:

pp = prey(0,[p_1;p_2])

% For the south-east isosector:

echo off;

p_1 = input('Pick a point on the horizontal axis in the isosector p_1=');
p_2 = input('Pick a point on the vertical axis in the isosector p_2=');

echo on;

   % The directions of the variables are:

pp = prey(0,[p_1;p_2])

% For the south-west isosector:

echo off;

p_1 = input('Pick a point on the horizontal axis in the isosector p_1=');
p_2 = input('Pick a point on the vertical axis in the isosector p_2=');

echo on;

   % The direction of the variables are:

pp = prey(0,[p_1;p_2])

% For the north-west isosector:

echo off;

p_1 = input('Pick a point on the horizontal axis in the isosector p_1=');
p_2 = input('Pick a point on the vertical axis in the isosector p_2=');

echo on;

   % The directions of the variables are:

pp = prey(0,[p_1;p_2])

% Finally, we can plot the system in phase space (y_1 and y_2);

plot(y(:,1),y(:,2));
title('Lotke-Volterra Equation - phase-plane plot');

% The graph of the phase plane (two dimensional) trajectories for both
% predator and prey depicts a closed orbit around the equilibrium point
% obtained by setting both differential equations equal to zero. The isoclines
% intersect at the equilibrium. By changing the parameter values and initial
% population size, it is possible to see different orbits.