% set parameters
a = 1.2;
b = 0.6;
c = 0.8;
d = 0.3;

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

% solve using Euler (eulsys.m)
[t, xeul] = eulsys(@predprey,[0 40],xini,h,a,b,c,d);

% solve using ODE45
[t, xode45] = ode45(@predprey,tspan,xini,[],a,b,c,d);

% plot
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')


% time step is too long for Euler
% try shorter time step
h = 0.00625;
[t2, xeul2] = eulsys(@predprey,[0 40],xini,h,a,b,c,d);

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)')