%
%  lotka.m	Run Gillespie example, eqn (38)
%
%  usage:	lotka(y,maxi,inc)
%
%		y = initial # of [y(1) y(2)] molecules
%		maxi = number of iterations to be taken
%		inc = number of iterations between displayed results
%
%  example:	lotka([2000 2000],100000,100)
%

function lotka(y,maxi,inc)

c1x = 10;
c2 = 0.01;
c3 = 10;

t = 0;

plot(y(1),y(2),'.')
hold on

for i=1:maxi,

    a1 = c1x*y(1);

    a2 = c2*y(1)*y(2);

    a3 = c3*y(2);

    a0 = a1 + a2 + a3;

    r1 = rand;
    
    tau = (1/a0)*log(1/r1);

    t = t + tau;

    r2 = rand;

    if r2 > (a1+a2)/a0

       	y(2) = y(2) - 1;	% apply reaction 3

    elseif r2 > a1/a0

        y(2) = y(2) + 1;	% apply reaction 2
        y(1) = y(1) - 1;

    else

        y(1) = y(1) + 1;	% apply reaction 1

    end

    if mod(i,inc) == 1

       plot(y(1),y(2),'.')

    end

end

hold off
xlabel('number of y_{1} molecules','fontsize',16)
ylabel('number of y_{2} molecules','fontsize',16)
title('Lotka via Gillespie','fontsize',16)