function Figure2_10a

% Matlab code for Figure 2.10(a) in book:
% Radek Erban and Jon Chapman, Stochastic Modelling of Reaction-Diffusion
% Processes, Cambridge University Press, 2019

close all;
rand('state',300);

k1=0.1;
k2=1;
k3=100;
k4=100;

for i=1:5

  A=0;
  B=20;
  C=0;
  time=0;
    
  k=1;
  while (k<500)
        rr=rand(2,1);
        a0=k1*A+k2*B*C+k3+k4*B*C;
        tau=(1/a0)*log(1/rr(1));
        time=time+tau;
        ss=k1*A;     
        if (ss>rr(2)*a0)
            A=A-1;
            B=B+1;
            k=k+1;
            Aplot(k,i)=A;
            timeplot(k,i)=time;
        else
            ss=ss+k2*B*C;
            if (ss>rr(2)*a0)
                A=A+1;
                B=B-1;
                k=k+1;
                Aplot(k,i)=A;
                timeplot(k,i)=time;
            else
                ss=ss+k3;
                if (ss>rr(2)*a0)
                   C=C+1;
                else
                   ss=ss+k4*B*C;
                   if (ss>rr(2)*a0)
                      C=C-1;
                   end
                end    
            end       
        end           
  end
end

[t,z] = ode15s(@myode,[0 2100],[0; 20; 0]);

t=t;

figure(1);
set(gca,'Fontsize',18);
plot(t,z(:,1),'--k','Linewidth',4);
hold on
h=stairs(timeplot(:,1),Aplot(:,1));
set(h,'Color','m','Linewidth',1);
h=stairs(timeplot(:,2),Aplot(:,2));
set(h,'Color','r','Linewidth',1);
h=stairs(timeplot(:,3),Aplot(:,3));
set(h,'Color','g','Linewidth',1);
h=stairs(timeplot(:,4),Aplot(:,4));
set(h,'Color','b','Linewidth',1);
h=stairs(timeplot(:,5),Aplot(:,5));
set(h,'Color','c','Linewidth',1);
plot(t,z(:,1),'--k','Linewidth',4);
xlabel('time [sec]','interpreter','latex');
ylabel('number of $A$ molecules','interpreter','latex');
h=legend('solution of ODEs');
set(h,'interpreter','latex','location','southeast','Fontsize',18);
axis([0 100 0 20]);
set(gca,'Fontsize',18);


function dydt = myode(t,z)
k1=0.1;
k2=1;
k3=10;
k4=10;
dydt =[k2*z(2)*z(3)-k1*z(1); k1*z(1)-k2*z(2)*z(3); k3-k4*z(2)*z(3)];

