% Predator-Prey-System
% Oliver Rheinbach 2010
% Universitaet Duisburg-Essen

function euler2
% Volterra Uebung 2.3
clear all;
close all;

% constants
alpha=2;
beta=1;
gamma=16;
%gamma=2;
%gamma=14;
%gamma=12;
delta=2;

% number of points
n=40000;

xl=0;
xr=4;
h=(xr-xl)/n;
x=xl:h:xr;
r=zeros(1,size(x,2));

% start values
if 1
b(1)=10;
r(1)=2;
else
b(1)=10;
r(1)=10;
end

for i=1:n
   b(i+1)=b(i) + h*(alpha*b(i) - beta*b(i)*r(i));
   r(i+1)=r(i) + h*(-gamma*r(i) + delta*b(i)*r(i));
end

% plot against the time
figure(1);
plot(x,b,'k-',x,r,'r-');

% plot r against b
ymin=min(r);
ymax=max(r);
xmin=min(b);
xmax=max(b);

figure(2);
vectorfield(xmin,xmax,ymin,ymax);

hold on
plot(b,r,'k');


function vectorfield(xmin,xmax,ymin,ymax);

% points
x=xmin:0.07*(xmax-xmin):xmax;
y=ymin:0.07*(ymax-ymin):ymax;

% make grid
[X,Y]=meshgrid(x,y);

F1=@(X,Y) [alpha*X - beta*(X.*Y)];
F2=@(X,Y) [-gamma*Y + delta*(Y.*X)];

U=F1(X,Y);
V=F2(X,Y);

% visualize
quiver(X,Y,U,V);
end

end