function Figure2_11a

% Matlab code for Figure 2.11(a) 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;

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

k5=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+k5*A*B+k5*A*B;
          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+k5*A*B;
                      if (ss>rr(2)*a0)
                         A=A+1;
                         B=B-1;                
                      else  
                         ss=ss+k5*A*B;
                         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;

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

k5=0.1;

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+k5*A*B+k5*A*B;
          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+k5*A*B;
                      if (ss>rr(2)*a0)
                         A=A+1;
                         B=B-1;                
                      else  
                         ss=ss+k5*A*B;
                         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;

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

k5=1;

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

p1=zeros(21,1);

for i=1:numberofrealisations
    
    while (time<i)
          rr=rand(2,1);
          a0=k1*A+k2*B*C+k3+k4*B*C+k5*A*B+k5*A*B;
          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+k5*A*B;
                      if (ss>rr(2)*a0)
                         A=A+1;
                         B=B-1;                
                      else  
                         ss=ss+k5*A*B;
                         if (ss>rr(2)*a0)
                            A=A-1;
                            B=B+1;
                         end
                      end
                   end
                end    
             end       
          end           
    end
    p1(A+1)=p1(A+1)+1;
end

p1=p1/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,p1);
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','Fontsize',18);
h=legend([h1 h2 h3],'$k_5=0$','$k_5 = 0.1$','$k_5 = 1$');
set(h,'interpreter','latex','location','north');
set(gca,'XTick',[0 2 4 6 8 10 12 14 16 18 20 22]);
axis([-0.5 20.5 0 1.05*max(p1)]);
set(gca,'Fontsize',18);



