% bimolecular decay of particles via Gillespie algorithm
%the reaction is S+S -> S
function stochastic_decay
global N alpha


set(0,                           ...
   'defaultaxesfontsize', 20,   ...
   'defaultaxeslinewidth', 1.0, ...
   'defaultlinelinewidth', 1.2, ...
   'defaultpatchlinewidth', 0.7);

%set the initial number of particles
N = 50; % initial particle number
% set the death rate alpha
alpha = 1;  % rate of decay
 
K = 10000; %total number of trials
 
t = zeros(N,K);% this is where we record the reaction times
R = rand(N,K); % generate the needed random numbers
for j = 1:N-1
    n = N-j+1;% number of particles left
    rxn = alpha*n*(n-1)/2; % reaction rate with n particles left
    dt = -log(R(j,:))/rxn; % next decay time interval
    t(j+1,:) = t(j,:) + dt; %update the time
    
    end
    s = 1; % use the first trial for plots
 tt = reshape([t(:,s)';t(:,s)'],2*(N),1);  % make a plot with stepwise decay
  nn = reshape([ N:-1:1; N:-1:1],2*(N),1);
%plot one decay trajectory
 figure(1)
 plot(tt(2:end),nn(1:end-1,1),tt,N./(1+N*alpha*tt/2),'--','linewidth',2)
 xlabel('t', 'fontsize', 20)
 ylabel('Particle Number', 'fontsize', 20)
 
 %plot one decay trajectory as algebraic decay 1/(1+a*t)
figure(2)
 plot(tt,1./nn,tt,(1/N+alpha*tt/2),'--','linewidth',2)
 xlabel('t', 'fontsize', 20)
 ylabel('1/Particle Number', 'fontsize', 20)
 title('1/Particle number vs. time','fontsize',20)
 
% create the pdf for extinction times from the data
  
 [NN,TT]=hist(t(N,:),50);  % histogram of the extinction times with 50 boxes
 tt = sort(TT);
   
 dt = mean(TT(2:end)-TT(1:end-1));  %timestep increment on the histogram
 
 % NN/(K*dt) is the approximate pdf for the extinction times from data
 % it is normalized to have total integral = 1
 
 %check the mean of the data
 mean(t(N,:))
 %vs the theoretical value
 2*(1-1/N)
 
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 % now simulate the differential equation system
 
 %dp_k/dt = a*(k+1)*p_(k+1) -a*k*p_k

% initialize the probabilities
p0 = [zeros(N-1,1);1]; %initialize p_k 
 tstep = 0.02;
t_end = 10;
tspan = [0:tstep:t_end];   
[T,S] = ode23s(@deRHS,tspan,p0); 
 
figure(3)
for j = 1:N
plot(T,S(:,j),'linewidth',2)
hold on
end
xlabel('t','fontsize',20)
ylabel('p_k(t)','fontsize',20)
hold off

% compare the data with the extinction time distribution dp_1/dt:
figure(4)
plot(TT,NN/(K*dt),'*',T,alpha*S(:,2),'--','linewidth',2)
 xlabel('\alpha t','fontsize',20)
ylabel('dp_1/dt','fontsize',20)

function p_prime=deRHS(t,p)  % right hand side for ode system
global N alpha 
 
K = [1:N]';
Km1 = [0:N-1]';
rt = K.*Km1/2;
p_prime = -alpha*rt.*p+alpha*[rt(2:end);0].*[p(2:end);0];

 

 
 
    