function prey_predator()
% =========================================================================
%  Biological_Models_Discrete_Prey_Predator.m
%  This is a script-function that simulates a discrete prey-predator model.
%  Note that an m-file can contain many functions.
% 
%  Pray - Predator Biological Models  
%
% This m-file describes a prey-predator systems.
% The general form of the system's equations is:
%
%  dN = 
%  dP = 
%
% Different models can be considered with different functional forms
%
% 04/15
% =========================================================================
 
% Clears history and command window 
close all; clear all;
 
% number of time steps
N = 150;
% -- const is a struct containing all constants for the model
const.K = 1;
%const.r = 1.3;
%const.s = 0.5;
%const.u = 0.7;
%const.v = 1.6;
const.r = 1.0;
const.s = 0.5;
const.u = 1.5;
const.v = 2.8;
 
% p: prey
% q: predator
p = zeros(1,N);
q = zeros(1,N);
 
% initial population
p(1) = 1.0; 
q(1) = 0.5;

% find populations at each time step
for i=2:N
    [p(i) q(i)] = F(p(i-1),q(i-1),const);  
end
 
% population vs time plot
figure (1);
%subplot(1,2,1)
plot(1:N,p(:),'r.-',1:N,q(:),'k.-')
legend('prey','predator')
xlabel('Time Steps')
ylabel('Population')
 
% phase plot
%subplot(1,2,2)
figure (2);
plot(p(:),q(:),'k.-')
xlabel('prey')
ylabel('predator')


%--------------------------------------------------------------------------
 
function [y1 y2] = F(x1,x2,c)
 
x1x2 = x1*x2;
 
y1 = x1 * (1+c.r*(1-x1/c.K)) - c.s * x1x2;
y2 =    (1-c.u)*x2           + c.v * x1x2;
