%%
clear
format compact

%%
alpha = 2/3;  
beta = 4/3;
gamma = 1;
delta = 1;
pp = @(t, u) predprey(t, u, alpha, beta, gamma, delta);

%%
% Note that the function must accept both t and u inputs, even though
% there is no explicit dependence on t. To solve the IVP we must also
% provide the initial condition, which is a 2-vector here, and the interval
% for the independent variable.
u0 = [1; 1];
[t, u] = eulersys(pp, [0, 40], u0, 6000);

%%
% Each row of the output u is the value of the solution vector at a
% single node in time. Each column of u is the entire history of one
% component of the solution.
size_u = size(u)
x = u(:,1);  
y = u(:,2);  
plot(t, x)
hold on
plot(t, y)
xlabel('t')
ylabel('u(t)')
title('Predator-prey solution')
legend('prey','predator')
grid on

%%
% When there are just two components, it's common to plot the solution in
% the _phase plane_, i.e., with u_1 and u_2 along the axes and time as
% a parameterization of the curve.
plot(x, y)
xlabel('x') 
ylabel('y')
title('Predator-prey phase plane')
grid on