function Figure7_1

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

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

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

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

A=25*ones(n,1);

f=0.001*[ones(n/2,1); -ones(n/2,1)];
kplus=D/(h*h)*ones(n,1)+f/(2*h);
kminus=D/(h*h)*ones(n,1)-f/(2*h);

kplus(n)=[];
kminus(1)=[];
time=0;
 
rr=rand(2,1);      
a0=dot(kplus,A(1:n-1))+dot(kminus,A(2:n));
time=time+(1/a0)*log(1/rr(1));

while (time<finaltime)
       ss=0;
       k=0;
       while ((ss<=rr(2)*a0)&(k<n-1))
              k=k+1;
              ss=ss+kplus(k)*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+kminus(k-1)*A(k);
           end
           A(k)=A(k)-1;
           A(k-1)=A(k-1)+1;
       end
       rr=rand(2,1);
       a0=dot(kplus,A(1:n-1))+dot(kminus,A(2:n));        
       time=time+(1/a0)*log(1/rr(1));
end

discretn=101;
dt=0.1;
histdet=25*ones(discretn,1);
fdet=0.001*[ones((discretn-1)/2,1); 0; -ones((discretn-1)/2,1)];
Sdet=0.00001*[0:(discretn-1)/2 (discretn-3)/2:-1:0]';
xdet=[0:(discretn-1)]/(discretn-1);
hh=dt*D*(discretn-1)*(discretn-1);
hhh=dt*(discretn-1)/2;
for t=0:dt:finaltime
    histdetold=histdet;
    ii=1;
       histdet(ii)=histdetold(ii)+hh*(histdetold(ii+1)-histdetold(ii))-2*hhh*(fdet(ii)*histdetold(ii));
    for ii=2:discretn-1
       histdet(ii)=histdetold(ii)+hh*(histdetold(ii+1)+histdetold(ii-1)-2*histdetold(ii))-hhh*(fdet(ii+1)*histdetold(ii+1)-fdet(ii-1)*histdetold(ii-1));
    end
    ii=discretn;
       histdet(ii)=histdetold(ii)+hh*(histdetold(ii-1)-histdetold(ii))+2*hhh*(fdet(ii)*histdetold(ii));
end

figure(1);
subplot(2,1,1);
set(gca,'Fontsize',18);
line(xdet,Sdet,'Color','r','Linewidth',3);
xlabel('$x$ [mm]','interpreter','latex');
ylabel('signal  $S(x)$','Color','r','interpreter','latex');
axis([0 1 0 0.00052]);
set(gca,'XTick',[0 0.2 0.4 0.6 0.8 1]);
set(gca,'YTick',[0 0.0002 0.0004]);
box on;
set(gca,'Fontsize',18);

subplot(2,1,2);
set(gca,'Fontsize',18);
line([0 0.5],[fdet(1) fdet(1)],'Color','b','Linewidth',3);
line([0.5 1],[fdet(discretn) fdet(discretn)],'Color','b','Linewidth',3);
xlabel('$x$ [mm]','interpreter','latex');
ylabel('$f(x)$  [mm$^{-1}$]','Color','b','interpreter','latex');
set(gca,'XTick',[0 0.2 0.4 0.6 0.8 1]);
axis([0 1 -0.0011 0.0011]);
box on;
set(gca,'Fontsize',18);

figure(2);
set(gca,'Fontsize',18);
hl1 = line([0:0.01:1],histdet,'Color','r','Linewidth',4);
axis([0 1 0 75]);
ax1 = gca;
xlabel('$x$ [mm]','interpreter','latex');
hold on;
hb=bar(mesh,A);
set(hb,'FaceColor',[0.88 0.88 0.88]);
ax1 = gca;
pos2=get(ax1,'Position');
pos2(1)=pos2(1)*1.04;
pos2(3)=pos2(3)*0.98;
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',4,'Parent',ax2);
set(ax1,'YLim',[0 75]);
set(ax1,'XTick',[0 0.2 0.4 0.6 0.8 1]);
set(ax1,'YTick',[0 10 20 30 40 50 60 70]);
set(ax2,'YLim',[0 75/25]);
set(ax2,'YTick',[0 1 2 3]);
set(gca,'Fontsize',18);
set(get(ax1,'YLabel'),'String','number per compartment','interpreter','latex');
set(get(ax1,'YLabel'),'Fontsize',18);
set(get(ax2,'YLabel'),'String','concentration [$\mu$m$^{-1}$]','interpreter','latex');
set(get(ax2,'YLabel'),'Fontsize',18);
hb2=legend('PDE (mean)');
set(hb2,'interpreter','latex','Fontsize',18,'Color',[1 1 1]);
set(get(ax1,'YLabel'),'Fontsize',18);
set(gca,'Fontsize',18);
