function Figure3_1d

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

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

k1=0.001;
k2=0.75;
k3=165;
k4=10000;
k5=200;

dt=0.01;
finaltime=800; 
n=finaltime/dt;

X=zeros(n+1,1);
time=[0:dt:finaltime];

for i=1:n
    X(i+1) = X(i)+dt*(-k1*X(i)*X(i)*X(i)+k2*X(i)*X(i)-k3*X(i)+k4)+k5*sqrt(dt)*randn(1,1);   
end

[t,z] = ode15s(@myode,[0 5],[0]);

figure(1);
set(gca,'Fontsize',18);
line([5000 5000],[0 0],'Color','b','Linewidth',4);
hold on;
line([5000 5000],[100 100],'Color','r','Linewidth',4);
plot(time,X,'b');   
hold on;
plot(t,z,'r','Linewidth',3);
line([5 800],[100 100],'Color','r','Linewidth',3);
xlabel('$t$','interpreter','latex');
ylabel('$X$','interpreter','latex');
set(get(gca,'ylabel'),'rotation',0);
hh=legend('SDE','ODE');
set(hh,'interpreter','latex','Fontsize',18,'location','northwest');
axis([0 800 0 500]);
box on;
set(gca,'Fontsize',18);

function dydt = myode(t,z)
k1=0.001;
k2=0.75;
k3=165;
k4=10000;
dydt = [-k1*z(1)*z(1)*z(1)+k2*z(1)*z(1)-k3*z(1)+k4];
