function Figure4_1

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

close all;
randn('state',10);

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

x=[-1:0.05:1];
y=[-1:0.05:1];
for i=1:41
    for j=1:41
        gauss(i,j)=(1/(4*D*pi*finaltime))*exp(-(x(i).*x(i)+y(j).*y(j))/(4*D*finaltime));
    end
end

figure(1);
set(gca,'Fontsize',17);
pcolor(x,y,gauss');

xlabel('$x$ [mm]','interpreter','latex');
ylabel('$y$ [mm]','interpreter','latex');
set(gca,'XTick',[-0.6 -0.4 -0.2 0 0.2 0.4 0.6]);
set(gca,'YTick',[-0.6 -0.4 -0.2 0 0.2 0.4 0.6]);
axis([-0.7 0.7 -0.7 0.7]);
shading flat;
view(0,90);
box on;
colormap jet;
colorbar;
shading interp;
set(gca,'Fontsize',17);

dt=10;
n=finaltime/dt;
numberofrealizations=1000000;
hist=zeros(40,40);

for i=1:numberofrealizations
    dX = sqrt(2*D*dt)*randn(1,n);   
    X = cumsum(dX);             
    dY = sqrt(2*D*dt)*randn(1,n);   
    Y = cumsum(dY);
    if (((X(n)>-1)&(X(n)<1))&((Y(n)>-1)&(Y(n)<1)))
            hist(round(20*X(n)+20.5),round(20*Y(n)+20.5))=hist(round(20*X(n)+20.5),round(20*Y(n)+20.5))+1;
    end
end

hist=hist/numberofrealizations/0.05/0.05;

xs=[-1.:0.05:0.95];
ys=[-1.:0.05:0.95];

figure(2);
set(gca,'Fontsize',17);
pcolor(xs,ys,hist');
xlabel('$x$ [mm]','interpreter','latex');
ylabel('$y$ [mm]','interpreter','latex');
set(gca,'XTick',[-0.6 -0.4 -0.2 0 0.2 0.4 0.6]);
set(gca,'YTick',[-0.6 -0.4 -0.2 0 0.2 0.4 0.6]);
axis([-0.7 0.7 -0.7 0.7]);
shading faceted;
view(0,90);
box on;
colormap jet;
colorbar;
set(gca,'Fontsize',17);

