function Figure4_6c

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

close all;
rand('state',100);

D=0.0001;          % 10^{-4} mm^2 sec^{-1}
finaltime=10*60;   % 10 minutes
numberofmolecules=2000;

n=40;
h=1/n;
mesh=[h/2:h:1-h/2];

A=ones(n,1)*numberofmolecules/40;

time=0;

while (time<finaltime)
       rr=rand(2,1);
       a0=2*D/(h*h)*(numberofmolecules-A(n)/2);          
       time=time+(1/a0)*log(1/rr(1));
       ss=0;
       k=0;
       while ((ss<=rr(2)*a0)&(k<n-1))
              k=k+1;
              ss=ss+D/(h*h)*A(k);
       end
       if (ss>rr(2)*a0) 
           A(k)=A(k)-1;
           A(k+1)=A(k+1)+1;
       else
           k=1;
           while ((ss<=rr(2)*a0)&(k<n))
                  k=k+1;
                  ss=ss+D/(h*h)*A(k);
           end
           if (ss>rr(2)*a0)
              A(k)=A(k)-1;
              A(k-1)=A(k-1)+1;
           else
              A(1)=A(1)-1;
              numberofmolecules=numberofmolecules-1;
           end
       end
end

boxes=101;
dt=0.1;
histdet=ones(boxes,1)*2000/40;
hh=dt*D*(boxes-1)*(boxes-1)/1/1;
for t=0:dt:finaltime
    histdetold=histdet;
    ii=1;
       histdet(ii)=0;
    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=10 min','Fontsize',18,'interpreter','latex');
xlabel('$x$ [mm]','interpreter','latex');
hold on;
h1=bar(mesh,A);
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,'PDE');
set(hh,'interpreter','latex','Fontsize',18,'location','southeast');
set(gca,'Fontsize',18);



