function Figure2_11b

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

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

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

%====================

k6=0;

A=10;
B=10;
C=0;
time=0;

p0=zeros(21,1);

for i=1:numberofrealisations
    
    while (time<i)
          rr=rand(2,1);
          a0=k1*A+k2*B*C+k3+k4*B*C+2*k6*A*(A-1)*(A-2)*(A-3)*(A-4)*(A-5)*(A-6)*(A-7)*(A-8)*(A-9)*B*(B-1)*(B-2)*(B-3)*(B-4)*(B-5)*(B-6)*(B-7)*(B-8)*(B-9);
          tau=(1/a0)*log(1/rr(1));
          time=time+tau;
          ss=k1*A;     
          if (ss>rr(2)*a0)
             A=A-1;
             B=B+1;
          else
             ss=ss+k2*B*C;
             if (ss>rr(2)*a0)
                A=A+1;
                B=B-1;
             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;                   
                   else
                      ss=ss+k6*A*(A-1)*(A-2)*(A-3)*(A-4)*(A-5)*(A-6)*(A-7)*(A-8)*(A-9)*B*(B-1)*(B-2)*(B-3)*(B-4)*(B-5)*(B-6)*(B-7)*(B-8)*(B-9);
                      if (ss>rr(2)*a0)
                         A=A+1;
                         B=B-1;                
                      else  
                         ss=ss+k6*A*(A-1)*(A-2)*(A-3)*(A-4)*(A-5)*(A-6)*(A-7)*(A-8)*(A-9)*B*(B-1)*(B-2)*(B-3)*(B-4)*(B-5)*(B-6)*(B-7)*(B-8)*(B-9);
                         if (ss>rr(2)*a0)
                            A=A-1;
                            B=B+1;
                         end
                      end
                   end
                end    
             end       
          end           
    end
    p0(A+1)=p0(A+1)+1;
end

p0=p0/numberofrealisations;

%====================

k6=0.0000000000001;

A=10;
B=10;
C=0;
time=0;

p01=zeros(21,1);

for i=1:numberofrealisations
    
    while (time<i)
          rr=rand(2,1);
          a0=k1*A+k2*B*C+k3+k4*B*C+2*k6*A*(A-1)*(A-2)*(A-3)*(A-4)*(A-5)*(A-6)*(A-7)*(A-8)*(A-9)*B*(B-1)*(B-2)*(B-3)*(B-4)*(B-5)*(B-6)*(B-7)*(B-8)*(B-9);
          tau=(1/a0)*log(1/rr(1));
          time=time+tau;
          ss=k1*A;     
          if (ss>rr(2)*a0)
             A=A-1;
             B=B+1;
          else
             ss=ss+k2*B*C;
             if (ss>rr(2)*a0)
                A=A+1;
                B=B-1;
             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;                   
                   else
                      ss=ss+k6*A*(A-1)*(A-2)*(A-3)*(A-4)*(A-5)*(A-6)*(A-7)*(A-8)*(A-9)*B*(B-1)*(B-2)*(B-3)*(B-4)*(B-5)*(B-6)*(B-7)*(B-8)*(B-9);
                      if (ss>rr(2)*a0)
                         A=A+1;
                         B=B-1;                
                      else  
                         ss=ss+k6*A*(A-1)*(A-2)*(A-3)*(A-4)*(A-5)*(A-6)*(A-7)*(A-8)*(A-9)*B*(B-1)*(B-2)*(B-3)*(B-4)*(B-5)*(B-6)*(B-7)*(B-8)*(B-9);
                         if (ss>rr(2)*a0)
                            A=A-1;
                            B=B+1;
                         end
                      end
                   end
                end    
             end       
          end           
    end
    p01(A+1)=p01(A+1)+1;
end

p01=p01/numberofrealisations;

%====================

k6=0.00000000001;

A=10;
B=10;
C=0;
time=0;

p02=zeros(21,1);

for i=1:numberofrealisations
    
    while (time<i)
          rr=rand(2,1);
          a0=k1*A+k2*B*C+k3+k4*B*C+2*k6*A*(A-1)*(A-2)*(A-3)*(A-4)*(A-5)*(A-6)*(A-7)*(A-8)*(A-9)*B*(B-1)*(B-2)*(B-3)*(B-4)*(B-5)*(B-6)*(B-7)*(B-8)*(B-9);
          tau=(1/a0)*log(1/rr(1));
          time=time+tau;
          ss=k1*A;     
          if (ss>rr(2)*a0)
             A=A-1;
             B=B+1;
          else
             ss=ss+k2*B*C;
             if (ss>rr(2)*a0)
                A=A+1;
                B=B-1;
             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;                   
                   else
                      ss=ss+k6*A*(A-1)*(A-2)*(A-3)*(A-4)*(A-5)*(A-6)*(A-7)*(A-8)*(A-9)*B*(B-1)*(B-2)*(B-3)*(B-4)*(B-5)*(B-6)*(B-7)*(B-8)*(B-9);
                      if (ss>rr(2)*a0)
                         A=A+1;
                         B=B-1;                
                      else  
                         ss=ss+k6*A*(A-1)*(A-2)*(A-3)*(A-4)*(A-5)*(A-6)*(A-7)*(A-8)*(A-9)*B*(B-1)*(B-2)*(B-3)*(B-4)*(B-5)*(B-6)*(B-7)*(B-8)*(B-9);
                         if (ss>rr(2)*a0)
                            A=A-1;
                            B=B+1;
                         end
                      end
                   end
                end    
             end       
          end           
    end
    p02(A+1)=p02(A+1)+1;
end

p02=p02/numberofrealisations;

nval=[0:1:20];
nval2=[-0.5:1:20.5];
figure(1);
set(gca,'Fontsize',18);
h1=plot(nval,p0,'ro','Linewidth',4,'Markersize',4);
hold on;
h2=stairs(nval2,[p01; p01(21)],'b','Linewidth',4);
h3=bar(nval,p02);
set(h3,'facecolor',[0.88 0.88 0.88]);
stairs(nval2,[p01; p01(21)],'b','Linewidth',4);
plot(nval,p0,'ro','Linewidth',4,'Markersize',4);
xlabel('number of molecules','interpreter','latex');
ylabel('stationary distribution','interpreter','latex');
h=legend([h1 h2 h3],'$k_6=0$','$k_6 = 10^{-13}$','$k_6 = 10^{-11}$');
set(h,'interpreter','latex','location','northeast','Fontsize',18);
set(gca,'XTick',[0 2 4 6 8 10 12 14 16 18 20 22]);
axis([-0.5 20.5 0 0.15]);
set(gca,'Fontsize',18);



