function Figure4_2b

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

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

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

while (time<finaltime)
       time=time+dt; 
       position=position+sqrt(2*D*dt)*randn(numberofmolecules,1);
       for i=1:numberofmolecules
           if (position(i)<0)
               position(i)=-position(i);
           end
           if (position(i)>=1)
               position(i)=2-position(i);
           end
       end
end

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

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

boxes=101;
histdet=zeros(boxes,1);
histdet(41)=numberofmolecules*100/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));
    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

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,56.5,'time=4 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(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 60]);
set(ax2,'YLim',[0 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,'mean');
set(hh,'interpreter','latex','Fontsize',18,'location','northeast');
set(gca,'Fontsize',18);


