%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Solve the Lotka-Volterra differential equations 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

clear all; close all;

k = [1, 0.1, 0.1];
LVdiffeq = @(t, Y) [k(1)*Y(1) - k(2)*Y(1)*Y(2);
                    k(2)*Y(1)*Y(2) - k(3)*Y(2)];

%% Time evolution plot

figure(1);

y0 = [4,10]';
[Tout, Yout] = ode45(LVdiffeq, [0 1000], y0);

plot(Tout, Yout, '-');
legend('X (Prey)','Y (Predator)');
xlabel('Time');
ylabel('Population Size');
axis([0,100,0,18])
grid on

%% Phase space plot

Told = Tout;
Yold = Yout;

figure(2);

for i = [1,2,4,6,10]
    y0 = [4, i]';
    [Tout, Yout] = ode45(LVdiffeq, [0 100], y0);

    figure(2);
    plot(Yout(:,1), Yout(:,2), '--'); 
    gtext(sprintf('(%d,%d)',y0(1),y0(2)))
    xlabel('X (Prey)');
    ylabel('Y (Predator)');
    hold on;
end

plot(Yold(:,1), Yold(:,2), 'LineWidth', 2); 