function Figure5_2

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

close all;
load data_Figure5_2_a.dat;
datatrajectory=reshape(data_Figure5_2_a,103,5,2);

Cplot(:,:)=datatrajectory(:,:,1);
timeplot(:,:)=datatrajectory(:,:,2);

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

load data_Figure5_1_d.dat
histA=data_Figure5_1_d(:,1);
histB=data_Figure5_1_d(:,2);
histC=data_Figure5_1_d(:,3);

figure(2);
set(gca,'Fontsize',18);
hb1=bar([0:59],histC);
set(hb1,'FaceColor',[0.88 0.88 0.88]);
axis([5.5 33.5 0 0.1]);
xlabel('$n$','interpreter','latex');
ylabel('probability that $C(50)=n$','interpreter','latex');

load data_Figure5_2_b.dat
histC2=data_Figure5_2_b(:,1);

hold on;
plot([0:59],histC2,'r','Linewidth',3);
hh1=legend('Gillespie','multiscale');
set(hh1,'interpreter','latex','Fontsize',18);
set(gca,'Fontsize',18);

load data_Figure5_2_c_d.dat;
histA2=data_Figure5_2_c_d(:,1);
histB2=data_Figure5_2_c_d(:,2);

figure(3);
set(gca,'Fontsize',18);
hb2=bar([0:59],histA);
set(hb2,'FaceColor',[0.88 0.88 0.88]);
axis([10.5 41.5 0 0.086]);
set(gca,'XTick',[15 20 25 30 35 40]);
xlabel('$n$','interpreter','latex');
ylabel('probability that $A(50)=n$','interpreter','latex');

hold on;
plot([0:59],histA2,'r','Linewidth',3);
hh2=legend('Gillespie','multiscale');
set(hh2,'interpreter','latex','Fontsize',18);
set(gca,'Fontsize',18);

figure(4);
set(gca,'Fontsize',18);
hb3=bar([0:59],histB);
set(hb3,'FaceColor',[0.88 0.88 0.88]);
axis([-0.5 20.5 0 0.17]);
xlabel('$n$','interpreter','latex');
ylabel('probability that $B(50)=n$','interpreter','latex');

hold on;
plot([0:59],histB2,'r','Linewidth',3);
hh3=legend('Gillespie SSA','multiscale SSA');
set(hh3,'interpreter','latex','Fontsize',18);
set(gca,'Fontsize',18);
