% Stochastic simulation of A+B<->C
% via Basic Stochastic Method

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

time = 2;
h=.001;

iterations = time/h;
t = 0;

for i=1:iterations
    p = rand;
    a1 = k1*a(i)*b(i)*h;
    a2 = k2*c(i)*h;
    if p < a1
        a(i+1) = a(i)-1;
        b(i+1) = b(i)-1; 
        c(i+1) = c(i)+1;  
    elseif p < a1+a2
        a(i+1) = a(i)+1;
        b(i+1) = b(i)+1; 
        c(i+1) = c(i)-1; 
    else
        a(i+1) = a(i);
        b(i+1) = b(i); 
        c(i+1) = c(i); 
    end 
    t(i+1) = t(i)+h;
end

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