function Figure7_3

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

close all;

[t,y]=ode45(@myode,[0 50],[0 1]);

figure(1);
subplot(2,1,1);
set(gca,'Fontsize',18);
line([0 5 5 30 30 50],[1 1 2 2 1.5 1.5],'Color','b','Linewidth',3);
xlabel('time [sec]','interpreter','latex');
ylabel('signal $C(t)$','interpreter','latex');
axis([0 50 0 2.5]);
box on;
set(gca,'Fontsize',18);

subplot(2,1,2);
set(gca,'Fontsize',18);
line(t,1+y(:,1),'Color','r','Linewidth',3);
xlabel('time [sec]','interpreter','latex');
ylabel('response $R(t)$','interpreter','latex');
axis([0 50 0 2.5]);
box on;
set(gca,'Fontsize',18);

D=0.0002;
h=0.01;
xdet=[-1+h/2:h:5-h/2];
w=size(xdet);
discretn=w(2);
Diff=D/(h*h);
detden=zeros(discretn,1);
detden(discretn/3)=1;
dt=0.1;
fdet=4*0.0004/3*[ones(2*(discretn)/3,1);  -ones((discretn)/3,1)];

for t=0:dt:60*60
    detdenold=detden;
    ii=1;     
       detden(ii)=detdenold(ii)+dt*Diff*(detdenold(ii+1)-detdenold(ii))-dt*(fdet(ii)*detdenold(ii))/(h);
    for ii=2:discretn-1
       detden(ii)=detdenold(ii)+dt*Diff*(detdenold(ii+1)+detdenold(ii-1)-2*detdenold(ii))-dt*(fdet(ii+1)*detdenold(ii+1)-fdet(ii-1)*detdenold(ii-1))/(2*h);
    end
    ii=discretn;
       detden(ii)=detdenold(ii)+dt*Diff*(detdenold(ii-1)-detdenold(ii))+dt*(fdet(ii)*detdenold(ii))/(h);
end

histogram=load('data_Figure7_3.dat');

figure(2);
set(gca,'Fontsize',18);
plot(xdet,4000*detden,'r','Linewidth',4);
hold on;
hb1=bar(histogram(:,2),histogram(:,1));
set(hb1,'FaceColor',[0.88 0.88 0.88]);
plot(xdet,4000*detden,'r','Linewidth',4);
box on;
xlabel('$x$ [mm]','interpreter','latex');
ylabel('number of bacteria','interpreter','latex');
hh1=legend('chemotaxis equation (7.36)');
set(hh1,'Fontsize',18,'interpreter','latex','location','northwest');
axis([0 4 0 50]);
set(gca,'Fontsize',18);

function dy=myode(t,y)

if (t<5)
    signal=1;
elseif (t<30)
    signal=2;
else
    signal=1.5;
end
dy=[(signal-y(1)-y(2))/0.01; signal-y(2)];


