function Figure4_2

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);
detden=zeros(boxdeterm,1);
influxrate=0.3;
k1=0.001;

for t=0:dt:10*60
    detdenold=detden;
    ii=1;     
       detden(ii)=detdenold(ii)+dt*Diff*(detdenold(ii+1)-detdenold(ii))-dt*k1*detdenold(ii);
    for ii=2:boxdeterm-1
       detden(ii)=detdenold(ii)+dt*Diff*(detdenold(ii+1)+detdenold(ii-1)-2*detdenold(ii))-dt*k1*detdenold(ii);
    end;    
    ii=boxdeterm;
       detden(ii)=detdenold(ii)+dt*Diff*(detdenold(ii-1)-detdenold(ii))-dt*k1*detdenold(ii);

    for ii=1:boxdeterm/5
       detden(ii)=detden(ii)+dt*influxrate;
    end;    
end;

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

figure(1);
set(gca,'Fontsize',20);
plot([0 xdet 1],[detden(1); detden; detden(boxdeterm)],'r','Linewidth',4);
hold on;
bar(mesh,histogram(:,1));
colormap([0.88 0.88 0.88]);
plot([0 xdet 1],[detden(1); detden; detden(boxdeterm)],'r','Linewidth',4);
box on;
xlabel('x [mm]');
ylabel('number of molecules');
text(0.5,95,'time=10 min','Fontsize',20);
legend('mean');
axis([0 1 0 100]);

for t=10*60+dt:dt:30*60
    detdenold=detden;
    ii=1;     
       detden(ii)=detdenold(ii)+dt*Diff*(detdenold(ii+1)-detdenold(ii))-dt*k1*detdenold(ii);
    for ii=2:boxdeterm-1
       detden(ii)=detdenold(ii)+dt*Diff*(detdenold(ii+1)+detdenold(ii-1)-2*detdenold(ii))-dt*k1*detdenold(ii);
    end;    
    ii=boxdeterm;
       detden(ii)=detdenold(ii)+dt*Diff*(detdenold(ii-1)-detdenold(ii))-dt*k1*detdenold(ii);

    for ii=1:boxdeterm/5
       detden(ii)=detden(ii)+dt*influxrate;
    end;    
end;

figure(2);
set(gca,'Fontsize',20);
plot([0 xdet 1],[detden(1); detden; detden(boxdeterm)],'r','Linewidth',4);
hold on;
bar(mesh,histogram(:,2));
colormap([0.88 0.88 0.88]);
plot([0 xdet 1],[detden(1); detden; detden(boxdeterm)],'r','Linewidth',4);
box on;
xlabel('x [mm]');
ylabel('number of molecules');
text(0.5,147.25,'time=30 min','Fontsize',20);
legend('mean');
set(gca,'YTick',[0 20 40 60 80 100 120 140]);
axis([0 1 0 155]);

