%% Lotka-Volterra Predator-prey Model
%
% This example demonstrates the use of different numerical ODE solvers 
% for a simple predator-pey model:
%    $$\frac{dx}{dt} = ax - bxy$
% and
%    $$\frac{dy}{dt} = -cy + dxy$
%
% First, set the model parameters:
a = 1.2;
b = 0.6;
c = 0.8;
d = 0.3;

%%
% Next, set initial values, time stepping etc.
h = 0.0625;
tspan = [0:h:40];
xini = [2 1];

%% Euler solution
%
% Solve using our own Euler solver (eulsys.m):
[t, xeul] = eulsys(@predprey,[0 40],xini,h,a,b,c,d);

%% Higher order Runge-Kutta solution
%
% Next solve using MATLAB's ODE45:
[t, xode45] = ode45(@predprey,tspan,xini,[],a,b,c,d);

%% Results
%
% Plot results:
figure
plot(t,xeul)
legend('prey','predator')
title('Euler solution')
xlabel('time')

figure
plot(t,xode45)
legend('prey','predator')
title('ODE45 solution')
xlabel('time')

figure
plot(xeul(1,:),xeul(2,:))
xlabel('prey')
ylabel('predator')
title('Euler solution')

figure
plot(xode45(:,1),xode45(:,2))
xlabel('prey')
ylabel('predator')
title('ODE45 solution')

%% Experimenting with different time steps
%
% It seems that our time step is too long for the Euler solver.
% Let's try shorter time step:
h = 0.00625;
[t2, xeul2] = eulsys(@predprey,[0 40],xini,h,a,b,c,d);

%%
% And plot...
figure
plot(t,xeul,t2,xeul2)
legend('prey (lg ts)','predator (lg ts)','prey (sm ts)','predator (sm ts)')

figure
plot(xeul2(1,:),xeul2(2,:))
xlabel('prey')
ylabel('predator')
title('Euler solution (small time step)')