% LV_gill.M
% Simple implementation of the Stochastic Simulation Algorithm
% (or Gillespie’s algorithm) for the Lotka-Volterra system.
%
% rand('state',100) % Can fix rand # pattern
% stoichiometric matrix for rxs
V = [1 -1 0;  0 1 -1];
%%%%%%% Parameters and Initial Conditions
X = zeros(2,1);
X(1) = 1000;   % initial molecules of X1
X(2) = 1000;   % initial molecules of X2
Y1(1) = X(1);  % store # of molecules
Y2(1) = X(2);
% set chem rx coefficients
c(1) = 10; c(2) = 0.01; c(3) = 10;
t = 0;         % initial time
T(1) = t;
tfinal = 10;   % final time
i = 1;
while t < tfinal
  % rx combination functions
  a(1) = c(1)*X(1);
  a(2) = c(2)*X(1)*X(2);
  a(3) = c(3)*X(2);
  asum = sum(a);  % total a
  % generate rand # and find rx occurring
  j = min(find(rand<cumsum(a/asum)));
  % 2nd rand # for time until rx
  tau = log(1/rand)/asum;
  X = X + V(:,j); % Stochastic matrix adjusts X
  t = t + tau;
  i = i + 1;
  T(i) = t;
  Y1(i) = X(1);
  Y2(i) = X(2);
end

figure(101)
plot(T,Y1,'b-',T,Y2,'r-');grid;
xlim([0,10]);
fontlabs = 'Times New Roman'; % Font type used in labels
xlabel('$t$','FontSize',14,'FontName',fontlabs,...
    'interpreter','latex');  
ylabel('Molecules','FontSize',14,'FontName',fontlabs);  
set(gca,'FontSize',12);

figure(102)
plot(Y1,Y2,'b-');grid;
xlabel('$Y_1$','FontSize',14,'FontName',fontlabs,...
    'interpreter','latex');  
ylabel('$Y_2$','FontSize',14,'FontName',fontlabs,...
    'interpreter','latex');
