function Figure7_2a

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

close all;

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

for t=0:dt:15*60
    detdenold=detden;
    ii=1;     
       detden(ii)=detdenold(ii)+dt*Diff*(detdenold(ii+1)-detdenold(ii))-dt*k1*detdenold(ii)-dt*(fdet(ii)*detdenold(ii))/(h);
    for ii=2:discretn-1
       detden(ii)=detdenold(ii)+dt*Diff*(detdenold(ii+1)+detdenold(ii-1)-2*detdenold(ii))-dt*k1*detdenold(ii)-dt*(fdet(ii+1)*detdenold(ii+1)-fdet(ii-1)*detdenold(ii-1))/(2*h);
    end
    ii=discretn;
       detden(ii)=detdenold(ii)+dt*Diff*(detdenold(ii-1)-detdenold(ii))-dt*k1*detdenold(ii)+dt*(fdet(ii)*detdenold(ii))/(h);

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

histogram=load('data_Figure7_2_a.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(discretn)],'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(discretn)],'k','Linewidth',4);
box on;
xlabel('$x$ [mm]','interpreter','latex');
ylabel('number of molecules','interpreter','latex');
set(gca,'XTick',[0 0.2 0.4 0.6 0.8 1]);
hh1=legend('PDE (mean)');
set(hh1,'Fontsize',18,'interpreter','latex');
axis([0 1 0 125]);
set(gca,'Fontsize',18);

