% Stochastic simulation of A+B<->C
% via Gillespei's Stochastic Algorithm
function gillespieABC


a(1) = 60;
b(1) = 50;
c(1) = 0;
k1 = .1;
k2 = 1;

time = 2;
t(1)=0;

i=1;
while ( t(i) < time )
    p = rand;
    a1 = k1*a(i)*b(i);
    a2 = k2*c(i);
    sa = a1 + a2;
    dt = -log(p)/sa;
    t(i+1)=t(i)+dt;
    p2 = rand;
    if p2 < a1/sa
        a(i+1) = a(i)-1;
        b(i+1) = b(i)-1; 
        c(i+1) = c(i)+1;  
    else
        a(i+1) = a(i)+1;
        b(i+1) = b(i)+1; 
        c(i+1) = c(i)-1; 
    end 
    i = i + 1;
end

plot(t,a); hold on;
plot(t,b,'r'); plot(t,c,'k');
legend('A','B','C')
xlabel('Time')
ylabel('Concentrations')