function Figure4_3

close all;

D=0.0001;
h=0.01;
xdet=[h/2:h:1-h/2];
w=size(xdet);
boxdeterm=w(2);
dt=0.01;
Diff=D/(h*h);
detdenA=zeros(boxdeterm,1);
detdenB=zeros(boxdeterm,1);
k1=0.001;
k2=0.01;
k3=1.2;
k4=1;


for t=0:dt:30*60
    detdenoldA=detdenA;
    detdenoldB=detdenB;
    ii=1;     
       detdenA(ii)=detdenoldA(ii)+dt*Diff*(detdenoldA(ii+1)-detdenoldA(ii))-2*dt*k1*detdenoldA(ii)*detdenoldA(ii)-dt*k2*detdenoldA(ii)*detdenoldB(ii);
       detdenB(ii)=detdenoldB(ii)+dt*Diff*(detdenoldB(ii+1)-detdenoldB(ii))-dt*k2*detdenoldA(ii)*detdenoldB(ii);
    for ii=2:boxdeterm-1
       detdenA(ii)=detdenoldA(ii)+dt*Diff*(detdenoldA(ii+1)+detdenoldA(ii-1)-2*detdenoldA(ii))-2*dt*k1*detdenoldA(ii)*detdenoldA(ii)-dt*k2*detdenoldA(ii)*detdenoldB(ii);
       detdenB(ii)=detdenoldB(ii)+dt*Diff*(detdenoldB(ii+1)+detdenoldB(ii-1)-2*detdenoldB(ii))-dt*k2*detdenoldA(ii)*detdenoldB(ii);
    end;    
    ii=boxdeterm;
       detdenA(ii)=detdenoldA(ii)+dt*Diff*(detdenoldA(ii-1)-detdenoldA(ii))-2*dt*k1*detdenoldA(ii)*detdenoldA(ii)-dt*k2*detdenoldA(ii)*detdenoldB(ii);
       detdenB(ii)=detdenoldB(ii)+dt*Diff*(detdenoldB(ii-1)-detdenoldB(ii))-dt*k2*detdenoldA(ii)*detdenoldB(ii);
    for ii=1:9*boxdeterm/10
       detdenA(ii)=detdenA(ii)+dt*k3;
    end;    
    for ii=2*boxdeterm/5:boxdeterm
       detdenB(ii)=detdenB(ii)+dt*k4;
    end;    
end;

histogram=load('dataFigure4_3.dat');
ww=size(histogram);
boxes=ww(1);
mesh=[1/(2*boxes):1/boxes:1-1/(2*boxes)]';

figure(1);
set(gca,'Fontsize',20);
h1=bar(mesh,histogram(:,1));
colormap([0.88 0.88 0.88]);
hold on;
h2=plot([0 xdet 1],[detdenA(1); detdenA; detdenA(boxdeterm)],'r','Linewidth',4);
box on;
xlabel('x [mm]');
ylabel('number of A molecules');
legend([h1 h2],'stochastic simulation','solution of PDEs');
axis([0 1 0 35]);

figure(2);
set(gca,'Fontsize',20);
h1=bar(mesh,histogram(:,2));
colormap([0.88 0.88 0.88]);
hold on;
h2=plot([0 xdet 1],[detdenB(1); detdenB; detdenB(boxdeterm)],'r','Linewidth',4);
box on;
xlabel('x [mm]');
ylabel('number of B molecules');
legend([h1 h2],'stochastic simulation','solution of PDEs',2);
axis([0 1 0 280]);