function Figure7_4

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

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

dt=0.01;
timesave=1;
eta=2;

N=10;
omega=randn(N,1);
time=0;
Uplot=zeros(1000,1);
timeplot=zeros(1000,1);

for i=1:1000
      while (time<i*timesave)
                 U=sum(omega)/N;
                 time=time+dt; 
                 omega=omega+(sign(U)-omega)*dt+eta*sqrt(dt)*randn(N,1);
      end
      Uplot(i)=sum(omega)/N;
      timeplot(i)=time;
end

figure(1);
set(gca,'Fontsize',18);
plot(timeplot,Uplot);
axis([0 1000*timesave -2 2]);
text(740*timesave,1.7,'$N=10$','Fontsize',18,'interpreter','latex');
xlabel('time $t$','interpreter','latex')
ylabel('average velocity $U(t)$','interpreter','latex');
set(gca,'Fontsize',18);

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

randn('state',5);

N=20;
omega=randn(N,1);
time=0;
Uplot2=zeros(1000,1);
timeplot2=zeros(1000,1);

for i=1:1000
      while (time<i*timesave)
                 U=sum(omega)/N;
                 time=time+dt; 
                 omega=omega+(sign(U)-omega)*dt+eta*sqrt(dt)*randn(N,1);
      end
      Uplot2(i)=sum(omega)/N;
      timeplot2(i)=time;
end

figure(2);
set(gca,'Fontsize',18);
plot(timeplot2,Uplot2);
axis([0 1000*timesave -2 2]);
text(800*timesave,1.7,'$N=20$','Fontsize',18,'interpreter','latex');
xlabel('time $t$','interpreter','latex')
ylabel('average velocity $U(t)$','interpreter','latex');
set(gca,'Fontsize',18);

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

N=10;
u=[-3:0.01:3];
statdist10=(sqrt(N/pi)/(eta*(1+erf(sqrt(N)/eta))))*exp((N/(eta*eta))*(-1+2*abs(u) - u.*u));
N=20;
statdist20=(sqrt(N/pi)/(eta*(1+erf(sqrt(N)/eta))))*exp((N/(eta*eta))*(-1+2*abs(u) - u.*u));
N=30;
statdist30=(sqrt(N/pi)/(eta*(1+erf(sqrt(N)/eta))))*exp((N/(eta*eta))*(-1+2*abs(u) - u.*u));

figure(3);
set(gca,'Fontsize',18);
plot(u,statdist10,'b','Linewidth',3);
hold on;
plot(u,statdist20,'r','Linewidth',3);
plot(u,statdist30,'k','Linewidth',3);
hh1=legend('$N=10$','$N=20$','$N=30$');
set(hh1,'Fontsize',18,'interpreter','latex','location','north');
xlabel('$u$','interpreter','latex')
ylabel('$p_s(u)$','rotation',0,'interpreter','latex');
axis([-2.3 2.3 0 0.8]);
set(gca,'Fontsize',18);

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

Nvec=[1:1:65];
integralform=(pi).*(1+erf(sqrt(Nvec)/eta)).*erfi(sqrt(Nvec)/eta);
asymptoticT=(2*eta*sqrt(pi)./(sqrt(Nvec))).*exp((Nvec/(eta*eta)));

figure(4);
set(gca,'Fontsize',18);
semilogy(Nvec,integralform,'b','Linewidth',3);
hold on;
semilogy(Nvec,asymptoticT,'r--','Linewidth',2);
xlabel('$N$','interpreter','latex')
ylabel('average switching time','interpreter','latex');
hh2=legend('equation (7.47)','equation (7.48)');
set(hh2,'Fontsize',18,'interpreter','latex','location','northwest');
set(gca,'Xtick',[10 20 30 40 50 60]);
set(gca,'Ytick',[1e1 1e2 1e3 1e4 1e5 1e6]);
axis([5 62 6 3*1e7/5]);
set(gca,'Fontsize',18);



