function Figure2_4

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

k1=0.001;
k2=0.01;
k3=1.2;
k4=1;

Ainitial=sqrt((k3-k4)/(2*k1));
Binitial=k4/(k2*Ainitial);
numberofrealisations=1000000;
pa=zeros(100,1);
p2D=zeros(30,30);

A=round(Ainitial);
B=round(Binitial);
time=0;

for i=1:numberofrealisations
   
    while (time<i)
          rr=rand(2,1);
          a0=k1*A*(A-1)+k2*A*B+k3+k4;
          time=time+(1/a0)*log(1/rr(1));
          ss=k1*A*(A-1);
          if (rr(2)*a0<ss)
              A=A-2;
          else
              ss=ss+k2*A*B;
              if (rr(2)*a0<ss)
                  A=A-1;
                  B=B-1;
              else
              ss=ss+k3;
                  if (rr(2)*a0<ss)
                      A=A+1;
                  else
                      B=B+1;
                  end;
              end;
          end;
    end;
    if ((A>0)&(A<101))
       pa(A)=pa(A)+1;
    end;
    if (((A>0)&(A<31))&((B>0)&(B<31)))
       p2D(A,B)=p2D(A,B)+1;
    end;
end;

pa=pa/numberofrealisations;
p2D=p2D/numberofrealisations;

figure(1);
set(gca,'Fontsize',20);
contourf(p2D',128);
xlabel('number of A molecules');
ylabel('number of B molecules');
axis([1 30 1 30 0 1.05*max(max(p2D))]);
shading flat;
view(0,90);
box on;
colorbar;

figure(2);
set(gca,'Fontsize',20);
bar(pa);
colormap([0.88 0.88 0.88]);
hold on;
bar([10],pa(10),'r');
xlabel('number of A molecules');
ylabel('stationary distribution');
set(gca,'XTick',[2 4 6 8 10 12 14 16 18 20 22 24]);
axis([0.5 25.5 0 1.05*max(pa)]);
