%
%  gill1.m	Run Gillespie example, eqn (29)
%
%  usage:	gill1(y,maxi,inc)
%
%		y = initial # of y molecules
%		maxi = number of iterations to be taken
%		inc = number of iterations between displayed results
%
%  example:	gill1(3000,30000,100)
%

function gill1(y,maxi,inc)

c1x = 5;
c2 = 0.005;

t = 0;

plot(t,y,'.')
hold on

for i=1:maxi,

    a1 = c1x*y;

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

    a0 = a1 + a2;

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

    t = t + tau;

    r2 = rand;

    if r2 > a1/a0

       	y = y - 2;	% apply reaction 2

    else

        y = y + 1;	% apply reaction 1

    end

    if mod(i,inc) == 1

       plot(t,y,'.')

    end

end

hold off
xlabel('time','fontsize',16)
ylabel('number of y molecules','fontsize',16)
title('Gillespie (29): x+y->2y->z, c_{1}x = 5, c_{2} = 0.005','fontsize',16)