function Figure4_7b

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

close all;
rand('state',5);
randn('state',5);

dt=0.1;
finaltime=10*60;   % 10 minutes
D=0.0001;          % 10^{-4} mm^2 sec^{-1
numberofmolecules=16000;
position=rand(numberofmolecules,1);
insolution=ones(numberofmolecules,1);
time=0;

K=0.001;
Pc=K*sqrt(pi*dt)/(2*sqrt(D));

while (time<finaltime)
       time=time+dt; 
       oldposition=position;
       position=position+sqrt(2*D*dt)*randn(numberofmolecules,1);
       for i=1:numberofmolecules
           if (insolution(i)>0.5)
              if (position(i)<0)
                  if (rand(1,1)<Pc)
                     insolution(i)=0;
                  else
                     position(i)=-position(i);
                  end   
              elseif (position(i)>=1)
                  position(i)=2-position(i);
              else
                  if (rand(1,1)<Pc*exp(-oldposition(i)*position(i)/(D*dt)))
                     insolution(i)=0;
                  end
              end
           end
       end
end

hist=zeros(40,1);
x=[0.0125:0.025:0.9875];

for i=1:numberofmolecules
    if (insolution(i)>0.5)
        hist(round(40*position(i)+0.5))=hist(round(40*position(i)+0.5))+1;
    end
end

boxes=101;
histdet=ones(boxes,1)*numberofmolecules/40;
hh=dt*D*(boxes-1)*(boxes-1)/1/1;
for t=0:dt:finaltime
    histdetold=histdet;
    ii=1;
       histdet(ii)=histdetold(ii)+hh*(histdetold(ii+1)-histdetold(ii))-dt*K*(boxes-1)*histdetold(ii);
    for ii=2:boxes-1
       histdet(ii)=histdetold(ii)+hh*(histdetold(ii+1)+histdetold(ii-1)-2*histdetold(ii));
    end
    ii=boxes;
       histdet(ii)=histdetold(ii)+hh*(histdetold(ii-1)-histdetold(ii));
end

iy=8;

figure(1);
set(gca,'Fontsize',18);
hl1 = line([0:0.01:1],histdet,'Color','r','Linewidth',3);
axis([0 1 0 60]);
ax1 = gca;
text(0.028,iy*56.5,'time=10 min','Fontsize',18,'interpreter','latex');
xlabel('$x$ [mm]','interpreter','latex');
hold on;
h1=bar(x,hist);
set(h1,'FaceColor',[0.88 0.88 0.88]);
pos2=get(ax1,'Position');
pos2(1)=pos2(1)*1.1;
pos2(3)=pos2(3)*0.95;
set(ax1,'Position',pos2);
ax2 = axes('Position',get(ax1,'Position'),'XTick',[],'XAxisLocation','top','YAxisLocation','right','Color','none','Fontsize',18);
hl2 = line([0:0.01:1],histdet/25,'Color','r','Linewidth',3,'Parent',ax2);
set(ax1,'YLim',[0 iy*60]);
set(ax2,'YLim',[0 iy*60/25]);
set(gca,'Fontsize',18);
set(get(ax1,'YLabel'),'String','molecules per compartment','interpreter','latex');
set(get(ax1,'YLabel'),'Fontsize',18);
set(get(ax2,'YLabel'),'String','concentration [molecules/$\mu$m]','interpreter','latex');
set(get(ax2,'YLabel'),'Fontsize',18);
hh=legend(hl1,'PDE');
set(hh,'interpreter','latex','Fontsize',18,'location','southeast');
set(gca,'Fontsize',18);

