function Figure6_1

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

close all;

D=0.0001;
h=0.01;
xdet=[h/2:h:1-h/2];
w=size(xdet);
boxdeterm=w(2);
dt=0.01;
Diff=D/(h*h);
detden=zeros(boxdeterm,1);
influxrate=0.3125;
k1=0.001;

for t=0:dt:10*60
    detdenold=detden;
    ii=1;     
       detden(ii)=detdenold(ii)+dt*Diff*(detdenold(ii+1)-detdenold(ii))-dt*k1*detdenold(ii);
    for ii=2:boxdeterm-1
       detden(ii)=detdenold(ii)+dt*Diff*(detdenold(ii+1)+detdenold(ii-1)-2*detdenold(ii))-dt*k1*detdenold(ii);
    end
    ii=boxdeterm;
       detden(ii)=detdenold(ii)+dt*Diff*(detdenold(ii-1)-detdenold(ii))-dt*k1*detdenold(ii);

    for ii=1:boxdeterm/5
       detden(ii)=detden(ii)+dt*influxrate;
    end
end

histogram=load('data_Figure6_1.dat');
ww=size(histogram);
boxes=ww(1);
mesh=[1/(2*boxes):1/boxes:1-1/(2*boxes)]';

figure(1);
set(gca,'Fontsize',18);
plot([0 xdet 1],[detden(1); detden; detden(boxdeterm)],'k','Linewidth',4);
hold on;
hb1=bar(mesh,histogram(:,1));
set(hb1,'FaceColor',[0.88 0.88 0.88]);
plot([0 xdet 1],[detden(1); detden; detden(boxdeterm)],'k','Linewidth',4);
box on;
xlabel('$x$ [mm]','interpreter','latex');
ylabel('number of molecules','interpreter','latex');
text(0.35,96.4,'time=10 min','Fontsize',18,'interpreter','latex');
set(gca,'XTick',[0 0.2 0.4 0.6 0.8 1]);
hh1=legend('mean');
set(hh1,'Fontsize',18,'interpreter','latex');
axis([0 1 0 105]);
set(gca,'Fontsize',18);

for t=10*60+dt:dt:30*60
    detdenold=detden;
    ii=1;     
       detden(ii)=detdenold(ii)+dt*Diff*(detdenold(ii+1)-detdenold(ii))-dt*k1*detdenold(ii);
    for ii=2:boxdeterm-1
       detden(ii)=detdenold(ii)+dt*Diff*(detdenold(ii+1)+detdenold(ii-1)-2*detdenold(ii))-dt*k1*detdenold(ii);
    end   
    ii=boxdeterm;
       detden(ii)=detdenold(ii)+dt*Diff*(detdenold(ii-1)-detdenold(ii))-dt*k1*detdenold(ii);

    for ii=1:boxdeterm/5
       detden(ii)=detden(ii)+dt*influxrate;
    end
end

figure(2);
set(gca,'Fontsize',18);
plot([0 xdet 1],[detden(1); detden; detden(boxdeterm)],'k','Linewidth',4);
hold on;
hb2=bar(mesh,histogram(:,2));
set(hb2,'FaceColor',[0.88 0.88 0.88]);
plot([0 xdet 1],[detden(1); detden; detden(boxdeterm)],'k','Linewidth',4);
box on;
xlabel('$x$ [mm]','interpreter','latex');
ylabel('number of molecules','interpreter','latex');
text(0.35,142,'time=30 min','Fontsize',18,'interpreter','latex');
hh2=legend('mean');
set(hh2,'Fontsize',18,'interpreter','latex');
set(gca,'XTick',[0 0.2 0.4 0.6 0.8 1]);
set(gca,'YTick',[0 20 40 60 80 100 120 140]);
axis([0 1 0 155]);
set(gca,'Fontsize',18);

