function Figure5_1

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

close all;
load data_Figure5_1_a_b_c.dat;
datatrajectory=reshape(data_Figure5_1_a_b_c,510,5,4);

Aplot(:,:)=datatrajectory(:,:,1);
Bplot(:,:)=datatrajectory(:,:,2);
Cplot(:,:)=datatrajectory(:,:,3);
timeplot(:,:)=datatrajectory(:,:,4);

figure(1);
set(gca,'Fontsize',18);
h=plot(timeplot(:,1),Aplot(:,1));
set(h,'Color','m','Linewidth',1);
hold on;
h=plot(timeplot(:,5),Aplot(:,5));
set(h,'Color','c','Linewidth',1);
h=plot(timeplot(:,2),Aplot(:,2));
set(h,'Color','g','Linewidth',1);
h=plot(timeplot(:,4),Aplot(:,4));
set(h,'Color','r','Linewidth',1);
h=plot(timeplot(:,3),Aplot(:,3));
set(h,'Color','b','Linewidth',1);
xlabel('time [sec]','interpreter','latex');
ylabel('number of $A$ molecules','interpreter','latex');
set(gca,'Fontsize',18);
axis([0 60 0 41]);

figure(2);
set(gca,'Fontsize',18);
h=plot(timeplot(:,1),Bplot(:,1));
set(h,'Color','m','Linewidth',1);
hold on;
h=plot(timeplot(:,5),Bplot(:,5));
set(h,'Color','c','Linewidth',1);
h=plot(timeplot(:,2),Bplot(:,2));
set(h,'Color','g','Linewidth',1);
h=plot(timeplot(:,4),Bplot(:,4));
set(h,'Color','r','Linewidth',1);
h=plot(timeplot(:,3),Bplot(:,3));
set(h,'Color','b','Linewidth',1);
xlabel('time [sec]','interpreter','latex');
ylabel('number of $B$ molecules','interpreter','latex');
set(gca,'Fontsize',18);
axis([0 60 0 16]);

figure(3);
set(gca,'Fontsize',18);
h=stairs(timeplot(:,1),Cplot(:,1));
set(h,'Color','m','Linewidth',1);
hold on;
h=stairs(timeplot(:,5),Cplot(:,5));
set(h,'Color','c','Linewidth',1);
h=stairs(timeplot(:,2),Cplot(:,2));
set(h,'Color','g','Linewidth',1);
h=stairs(timeplot(:,4),Cplot(:,4));
set(h,'Color','r','Linewidth',1);
h=stairs(timeplot(:,3),Cplot(:,3));
set(h,'Color','b','Linewidth',1);
xlabel('time [sec]','interpreter','latex');
ylabel('number of $C$ molecules','interpreter','latex');
set(gca,'Fontsize',18);
axis([0 60 0 33]);

load data_Figure5_1_d.dat;
histC=data_Figure5_1_d(:,3);
[t,y]=ode45(@detODE,[0 50],[0; 0; 0]);
lambda=y(end,3);
nvalues=[5:34];
poissonsol=exp(-lambda)*(lambda.^nvalues)./factorial(nvalues);

figure(4);
set(gca,'Fontsize',18);
hb=bar([0:59],histC);
set(hb,'FaceColor',[0.88 0.88 0.88]);
hold on;
plot(nvalues,poissonsol,'ro','Linewidth',4,'Markersize',3);
hh1=legend('Gillespie','Poisson');
set(hh1,'interpreter','latex','Fontsize',18);
axis([5.5 33.5 0 0.1]);
xlabel('$n$','interpreter','latex');
ylabel('probability that $C(50)=n$','interpreter','latex');
set(gca,'Fontsize',18);

function dydt=detODE(t,y)
dydt=[-50 200 0; 50 -200.1 0; 0 0.1 0]*y+[1; 0; 0];
