function Figure3_3a

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

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

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

A=zeros(n,1);

A(15)=numberofmolecules/2;
A(16)=numberofmolecules/2;

time=0;

while (time<finaltime)
       rr=rand(2,1);
       a0=2*D/(h*h)*(numberofmolecules-A(1)/2-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;
           A(k)=A(k)-1;
           A(k-1)=A(k-1)+1;
       end;
end;

boxes=101;
dt=0.1;
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',20);
hl1 = line([0:0.01:1],histdet,'Color','r','Linewidth',4);
ax1 = gca;
text(0.028,56.5,'time=4 min','Fontsize',20);
legend('mean');
xlabel('x [mm]');
hold on;
bar(mesh,A);
colormap([0.88 0.88 0.88]);
ax2 = axes('Position',get(ax1,'Position'),'XTick',[],'XAxisLocation','top','YAxisLocation','right','Color','none','Fontsize',20);
hl2 = line([0:0.01:1],histdet/25,'Color','r','Linewidth',4,'Parent',ax2);
set(ax1,'YLim',[0 60]);
set(ax2,'YLim',[0 60/25]);
set(get(ax1,'YLabel'),'String','number of molecules in compartment');
set(get(ax1,'YLabel'),'Fontsize',20);
set(get(ax2,'YLabel'),'String','concentration in [molecules/\mu{}m]');
set(get(ax2,'YLabel'),'Fontsize',20);


