function Figure2_1

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

close all;

[t1,z1] = ode45(@myode,[0 1],[0]);
[t2,z2] = ode45(@myode,[0 1],[200]);
[t3,z3] = ode45(@myode,[0 1],[300]);
[t4,z4] = ode45(@myode,[0 1],[500]);

figure(1);
set(gca,'Fontsize',18);
plot(t1,z1,'r','Linewidth',3);
hold on;
plot(t2,z2,'k','Linewidth',3);
plot(t3,z3,'g','Linewidth',3);
plot(t4,z4,'b','Linewidth',3);
line([0 1],[100 100],'LineStyle',':','Color','k');
hold on;
line([0 1],[220 220],'LineStyle',':','Color','k');
line([0 1],[400 400],'LineStyle',':','Color','k');
xlabel('time [min]','interpreter','latex');
ylabel('number of $A$ molecules','interpreter','latex');
hh=legend('$A(0)=0$','$A(0)=200$','$A(0)=300$','$A(0)=500$');
set(hh,'interpreter','latex','location','east','Fontsize',18);
axis([0 1 0 550]);
box on;
set(gca,'Fontsize',18);

rand('state',44);

k1=0.00025;
k2=0.18;
k3=37.5;
k4=2200;

X=0;
time=0;
timeSSA=0;
kk=0;

while (time<2)
   timeSSAprev=timeSSA; 
   timefinSSA=1;
   kk=kk+1;
   while (timeSSA<timeSSAprev+timefinSSA)
      rr=rand(2,1);
      timeSSA=timeSSA+1;
      a0=k1*X*(X-1)*(X-2)+k2*X*(X-1)+k3*X+k4;
      tau=(1/a0)*log(1/rr(1));
      if (rr(2)*a0<(k2*X*(X-1)+k4))
          X=X+1;
      else
          X=X-1;
      end
      time=time+tau;
   end;     
XX(kk)=X;
tt(kk)=time;
end

[t,z] = ode45(@myode,[0 2],[0]);

figure(2);
set(gca,'Fontsize',18);
line([5 5],[0 0],'Color','b','Linewidth',4);
hold on;
line([5 5],[0 0],'Color','r','Linewidth',4);
h=stairs(tt,XX);
set(h,'Color','b','Linewidth',1);
plot(t,z,'r','Linewidth',3);
xlabel('time [min]','interpreter','latex');
ylabel('number of $A$ molecules','interpreter','latex');
hh=legend('stochastic','deterministic');
set(hh,'interpreter','latex','Fontsize',18);
axis([0 2 0 550]);
box on;
set(gca,'Fontsize',18);

function dydt = myode(t,z)
k1=0.00025;
k2=0.18;
k3=37.5;
k4=2200;
dydt = [-k1*z(1)*z(1)*z(1)+k2*z(1)*z(1)-k3*z(1)+k4];
