function Figure3_4a

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

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

k1=0.1;
k2=1;

A=k2/k1;
numberofrealisations=100000;
time=0;

p=zeros(22,1);

for i=1:numberofrealisations
    
    while (time<i)
          rr=rand(2,1);
          a0=k1*A+k2;
          time=time+(1/a0)*log(1/rr(1));
          if (rr(2)*a0<k1*A)
              A=A-1;
          else
              A=A+1;
          end
    end
    if ((A>0)&(A<23))
       p(A)=p(A)+1;
    end
end

p=p/numberofrealisations;

av=[0:1:25];
pFP=exp(-2*av+((4*k2/k1)-1)*log(k1*av+k2));
pFP=pFP/sum(pFP);

figure(1);
set(gca,'Fontsize',17);
h1=bar(p);
set(h1,'FaceColor',[0.88 0.88 0.88]);
hold on;
h2=plot(av,pFP,'r','Linewidth',3);
xlabel('number of molecules','interpreter','latex');
ylabel('stationary distribution','interpreter','latex');
h3=legend([h1 h2 ],'$\;$Gillespie SSA','$\;$Fokker-Planck');
set(h3,'interpreter','latex','Fontsize',17,'location','northeast');
set(gca,'XTick',[0 2 4 6 8 10 12 14 16 18 20 22]);
axis([0 22.5 0 1.05*max(p)]);
set(gca,'Fontsize',17);

