function Figure2_2b

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;

pdet(1)=1;
pdet(2)=pdet(1)*1/(2)+(k2/k1)*(1/(2))*(pdet(1));
int=pdet(1)+pdet(2);
for i=2:25
pdet(i+1)=pdet(i)*i/(i+1)+(k2/k1)*(1/(i+1))*(pdet(i)-pdet(i-1));
int=int+pdet(i+1);
end;
pdet=pdet/int;

figure(1);
set(gca,'Fontsize',20);
h1=bar(p);
colormap([0.88 0.88 0.88]);
hold on;
h2=plot(pdet,'r','Linewidth',4);
xlabel('number of molecules');
ylabel('stationary distribution');
legend([h1 h2],'Gillespie SSA','master equation');
set(gca,'XTick',[0 2 4 6 8 10 12 14 16 18 20 22]);
axis([0 22.5 0 1.05*max(p)]);