function Figure2_8and2_9b

% Matlab code for Figures 2.8 and 2.9(b) in book:
% Radek Erban and Jon Chapman, Stochastic Modelling of Reaction-Diffusion
% Processes, Cambridge University Press, 2019

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

k1=100;
k2=0.1;    
k3=0.01;  
k4=k2*99/10;   
k5=1000;
k6=100;

for i=1:5

  A=10;
  B=100;
  C=10;
  time=0;
  k5=1000;
    
  for kk=1:2100
      while (time<kk*1)
           if (time>600) 
               k5=500;
           end 
           rr=rand(2,1);
           a0=k1+k2*C+k3*B+k4*A*C+k5+k6*A;
           tau=(1/a0)*log(1/rr(1));
           ss=k1;     
           if (ss>rr(2)*a0)
              C=C+1;
           else
              ss=ss+k2*C;
              if (ss>rr(2)*a0)
                 C=C-1;
                 B=B+1;
              else
                 ss=ss+k3*B;
                 if (ss>rr(2)*a0)
                    B=B-1;
                 else
                    ss=ss+k4*A*C;
                    if (ss>rr(2)*a0)
                       C=C-1;
                    else
                       ss=ss+k5;
                       if (ss>rr(2)*a0)
                          A=A+1;
                       else
                          ss=ss+k6*A; 
                          if (ss>rr(2)*a0)
                             A=A-1;
                          end
                       end
                    end
                 end
              end
           end
           time=time+tau;
      end
      AA(kk,i)=A;
      BB(kk,i)=B;
      CC(kk,i)=C;
      tt(kk,i)=time/60;
  end
end

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

t=t/60;

figure(1);
set(gca,'Fontsize',18);
plot(t,z(:,1),'k','Linewidth',4);
hold on;
plot(tt(:,1),AA(:,1),'b');
plot(tt(:,1),AA(:,2),'r');
plot(tt(:,1),AA(:,3),'g');
plot(tt(:,1),AA(:,4),'m');
plot(tt(:,1),AA(:,5),'y');
plot(t,z(:,1),'k','Linewidth',4);
xlabel('time [min]','interpreter','latex');
ylabel('number of $A$ molecules','interpreter','latex');
hh=legend('solution of ODEs');
set(hh,'interpreter','latex','location','northeast','Fontsize',18);
axis([0 35 0 22]);
set(gca,'Fontsize',18);

figure(2);
set(gca,'Fontsize',18);
plot(t,z(:,2),'k','Linewidth',4);
hold on;
plot(tt(:,1),BB(:,1),'b');
plot(tt(:,1),BB(:,2),'r');
plot(tt(:,1),BB(:,3),'g');
plot(tt(:,1),BB(:,4),'m');
plot(tt(:,1),BB(:,5),'y');
plot(t,z(:,2),'k','Linewidth',4);
xlabel('time [min]','interpreter','latex');
ylabel('number of $B$ molecules','interpreter','latex');
hh=legend('solution of ODEs');
set(hh,'interpreter','latex','location','northwest','Fontsize',18);
axis([0 35 50 400]);
set(gca,'Fontsize',18);

figure(3);
set(gca,'Fontsize',18);
plot(t,z(:,3),'k','Linewidth',4);
hold on;
plot(tt(:,1),CC(:,1),'b');
plot(tt(:,1),CC(:,2),'r');
plot(tt(:,1),CC(:,3),'g');
plot(tt(:,1),CC(:,4),'m');
plot(tt(:,1),CC(:,5),'y');
plot(t,z(:,3),'k','Linewidth',4);
xlabel('time [min]','interpreter','latex');
ylabel('number of $C$ molecules','interpreter','latex');
hh=legend('solution of ODEs');
set(hh,'interpreter','latex','location','southeast','Fontsize',18);
axis([0 35 0 max(max(CC))+1]);
set(gca,'Fontsize',18);

function dydt = myode(t,z)
k1=100;
k2=0.1;
k3=0.01;
k4=k2*99/10;
k5=1000;
k6=100;
if (t>600)
   k5=500;
end   
dydt =[k5-k6*z(1); k2*z(3)-k3*z(2); k1-k2*z(3)-k4*z(1)*z(3)];

