function Figure4_4_a_b

% Matlab code for Figures 4.4(a) and 4.4(b) in book:
% Radek Erban and Jon Chapman, Stochastic Modelling of Reaction-Diffusion
% Processes, Cambridge University Press, 2019

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

dt=0.00000000000002;
finaltime=0.0000000006;  % 10 minutes
D=0.0001;         % 10^{-4} mm^2 sec^{-1}
speed=10000;
lambda=speed*speed/(2*D);

n=finaltime/dt+1;
numberofrealisations=6;
X=zeros(numberofrealisations,n);

for i=1:numberofrealisations
    X(i,1)=0;
    if (rand(1,1)<0.5)
       velocity=speed;
    else
       velocity=-speed;
    end  
    for j=2:n
        if (rand(1,1)<lambda*dt)
            velocity=-velocity;
        end
        X(i,j)=X(i,j-1)+velocity*dt;
    end
end

time=[0:dt:finaltime]*1000*1000*1000;
X=X*1000*1000;

figure(1);
set(gca,'Fontsize',18);
plot(time,X(1,:),'r');
hold on;
plot(time,X(2,:),'c');
plot(time,X(3,:),'g');
plot(time,X(4,:),'b');
plot(time,X(5,:),'m');
plot(time,X(6,:),'Color',[0.9 0.9 0]);
xlabel('time [ns]','interpreter','latex');
ylabel('$x$ [nm]','interpreter','latex');
set(gca,'XTick',[0 0.1 0.2 0.3 0.4 0.5 0.6]);
set(gca,'YTick',[-0.5 -0.4 -0.3 -0.2 -0.1 0 0.1 0.2 0.3 0.4 0.5]);
axis([0 0.6 -0.5 0.5]);
set(gca,'Fontsize',18);

%=====================================

rand('state',80);

D=0.0001;         % 10^{-4} mm^2 sec^{-1}
speed=10000;
lambda=speed*speed/(2*D);

n=334;
X2=zeros(numberofrealisations,n);
time2=zeros(numberofrealisations,n);

for i=1:numberofrealisations
    if (rand(1,1)<0.5)
       velocity=speed;
    else
       velocity=-speed;
    end
    for j=2:n
        tau=(1/lambda)*log(1/rand(1,1));
        time2(i,j)=time2(i,j-1)+tau;
        X2(i,j)=X2(i,j-1)+velocity*tau;
        velocity=-velocity;
    end
end

time2=time2*1000*1000*1000;
X2=X2*1000*1000;

figure(2);
set(gca,'Fontsize',18);
plot(time2(1,:),X2(1,:),'r');
hold on;
plot(time2(2,:),X2(2,:),'b');
plot(time2(3,:),X2(3,:),'c');
plot(time2(4,:),X2(4,:),'m');
plot(time2(5,:),X2(5,:),'g');
plot(time2(6,:),X2(6,:),'Color',[0.9 0.9 0]);
xlabel('time [ns]','interpreter','latex');
ylabel('$x$ [nm]','interpreter','latex');
set(gca,'XTick',[0 0.1 0.2 0.3 0.4 0.5 0.6]);
set(gca,'YTick',[-0.5 -0.4 -0.3 -0.2 -0.1 0 0.1 0.2 0.3 0.4 0.5]);
axis([0 0.6 -0.5 0.5]);
set(gca,'Fontsize',18);
