function Figure6_6

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

close all;

% ===========================
% Figure 6.6(a)

load data_Figure6_6_a.dat
kammagamma05=data_Figure6_6_a(:,1);
kammagamma1=data_Figure6_6_a(:,2);
kammagamma2=data_Figure6_6_a(:,3);
Plambda=data_Figure6_6_a(:,4);

figure(1);
set(gca,'Fontsize',18);
g3=plot(kammagamma2,Plambda,'k--','Linewidth',3);
hold on;
g2=plot(kammagamma1,Plambda,'k','Linewidth',3);
g1=plot(kammagamma05,Plambda,'k:','Linewidth',3);
%plot([0:0.1:4.2],[0:0.1:4.2]*3/(4*pi),'g--');
hl=legend([g1 g2 g3],'$\gamma = 0.5$','$\gamma = 1$','$\gamma = 2$');
set(hl,'interpreter','latex','Fontsize',18,'location','southeast');
axis([0 3.9 0 1]);
xlabel('$\kappa$','interpreter','latex');
text(-0.6,0.9,'$P_\lambda$','Fontsize',18,'interpreter','latex');
set(gca,'XTick',[0 0.5 1 1.5 2 2.5 3 3.5]);
set(gca,'Fontsize',18);

% ===========================
% Figure 6.6(b)

K=0.002;
DiffConst=0.0001;
Kd=0.005;

h=0.02;
xdet=[0:h:1];
w=size(xdet);
boxdeterm=w(2);
dt=0.1;
Diff=DiffConst/(h*h);
detden=[ones((boxdeterm-1)/2,1); 0.5; zeros((boxdeterm-1)/2,1)];
bounded=0;

for t=0:dt:600
    boundedold=bounded; 
    detdenold=detden;
    ii=1;     
       bounded=boundedold+dt*(K*detdenold(ii)-Kd*boundedold);
       detden(ii)=1/(1+K*h/DiffConst)*detdenold(ii+1)+Kd*h/DiffConst/(1+K*h/DiffConst)*boundedold;
    for ii=2:boxdeterm-1
       detden(ii)=detdenold(ii)+dt*Diff*(detdenold(ii+1)+detdenold(ii-1)-2*detdenold(ii));
    end
    ii=boxdeterm;
       detden(ii)=detdenold(ii)+dt*Diff*(2*detdenold(ii-1)-2*detdenold(ii));
end

detden=detden*20000/40;

histogram=load('data_Figure6_6_b.dat');
ww=size(histogram);
boxes=ww(1);

figure(2);
set(gca,'Fontsize',18);
hb1=bar([1/(2*boxes):1/boxes:1-1/(2*boxes)]',[histogram]/25);
set(hb1,'FaceColor',[0.88 0.88 0.88]);
xlabel('$x$ [mm]','interpreter','latex');
ylabel('concentration [$\mu$m$^{-1}$]','interpreter','latex');
hold on;
colormap([0.88 0.88 0.88]);
hl=plot(xdet,detden/25,'k','Linewidth',4);
hh2=legend(hl,'PDE');
set(hh2,'interpreter','latex','Fontsize',18);
axis([0 1 0 270/25]);
set(gca,'Fontsize',18);

         