% Stochastic simulation of A+B<->C
% via Chemical Langevin Equation

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

time = 2;
dt=.01;

iterations = time/dt;

V = [ -1 1; -1 1; 1 -1 ];

t = 0;
for i=1:iterations
    prop = [ k1*a(i)*b(i); k2*c(i) ];
    x = [ a(i); b(i); c(i) ];
    y = x + V * prop * dt + V * ( sqrt(prop) .* randn(2,1) ) .* sqrt(dt);
    a(i+1) = y(1);
    b(i+1) = y(2);
    c(i+1) = y(3);
    t(i+1) = t(i) + dt;
end

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