function prey_predator()
% This is a script-function that simulates the prey-predator model.
% Note that an m-file can contain many functions.

clc; clear;

% number of time steps
N = 300;
% const is a struct containing all constants for the model
const.K = 1;
const.r = 1;
const.s = 0.5;
const.u = 0.7;
const.v = 1.6;

% p: prey
% q: predator
p = zeros(1,N);
q = zeros(1,N);

% initial population
p(1) = 1.1; 
q(1) = 0.4;

for i=2:N
    [p(i) q(i)] = F(p(i-1),q(i-1),const);  
end

% population vs time plot
subplot(1,2,1)
plot(1:N,p(:),'r.-',1:N,q(:),'k.-')
legend('prey','predator')

% phase plot
subplot(1,2,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;

