%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Stochastically simulate a Lotka-Volterra system
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

clear all; close all;

%% Define the reaction

% Y1 -> 2Y1
% Y1 + Y2 -> 2Y2
% Y2 -> 0

lhs = [1, 0;
       1, 1;
       0, 1];
   
rhs = [2, 0;
       0, 2;
       0, 0];
   
c = [1, 0.005, 0.6]';
y00 = [50, 100]';

T_MAX = 30;

%% Run the simulation

q = 1;
for y0 = [2.^(3:8);2*(2.^(3:8))]
    
sim_log = ssas_gillespie(lhs, rhs, c, y0, T_MAX, 0.01);

%%% Time evolution plot

figure(1);
subplot(2,3,q);
plot(sim_log(:,1), sim_log(:,2:end), '-');
legend('X (Prey)','Y (Predator)');
xlabel('Time');
xlim([0 30]);
ylim([0 2000]);
title(sprintf('X(0) = %d, Y(0) = %d', y0));
ylabel('Population Size');
q = q + 1
end

%% Phase space plot

figure(2);

plot(sim_log(:,2), sim_log(:,3), '-'); 
xlabel('X (Prey)');
ylabel('Y (Predator)');