% exponential decay of particles via Gillespie algorithm
%Reaction
%S-> 0 at rate alpha
% Use this code as a template for death/decay processes
function stochastic_decay
global N alpha

set(0,                           ...
   'defaultaxesfontsize', 20,   ...
   'defaultaxeslinewidth', 1.2, ...
   'defaultlinelinewidth', 2.0, ...
   'defaultpatchlinewidth', 0.7);
 
%set the initial number of particles
N = 25; % initial particle number
% set the death rate alpha
alpha = 1;  % rate of decay
 
K = 10000; %total number of trials
 
t = zeros(N+1,K);% this is where we record the reaction times
R = rand(N,K); % generate the needed random numbers.  There are exactly N steps for K trials
for j = 1:N
    n = N-j+1;% number of particles remaining
    rxn = alpha*n; % 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
    js = 1; % use the first trial for plots
 tt = reshape([t(:,js)';t(:,js)'],2*(N+1),1);  % make a plot with stepwise decay
  nn = reshape([ N:-1:0; N:-1:0],2*(N+1),1);
%plot one decay trajectory
 figure(1)
 plot(tt(2:end),nn(1:end-1,1),tt,N*exp(-tt),'--','linewidth',2)
 xlabel('\alpha t', 'fontsize', 20)
 ylabel('Particle Number', 'fontsize', 20)
  
 
 %plot one decay trajectory on a log scale
figure(2)
 semilogy(tt,nn,tt,N*exp(-tt),'--','linewidth',2)
 xlabel('\alpha t', 'fontsize', 20)
 ylabel('Particle Number', 'fontsize', 20)
 title('semilog plot: Particle number vs. time','fontsize',20)
 
% create the pdf for extinction times from the data
  
 [NN,TT]=hist(t(N+1,:),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
 
 % the sample mean:
 Sample_mean=mean(t(N+1,:))
 % compare with the theoretical mean:
 E = sum(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]; %initialize p_k % p_N(0) = 1
 tstep = 0.02;
t_end = 12;
tspan = [0:tstep:t_end];   
[T,S] = ode23s(@deRHS,tspan,p0); 
 
figure(3)  % plot the solution of p_k
for j = 1:N+1
plot(T,S(:,j),'linewidth',2)
hold on
end
xlabel('\alpha t','fontsize',20)
ylabel('p_k(t)','fontsize',20)
hold off

% compare the data with the extinction time distribution dp_0/dt:
figure(4)
plot(TT,NN/(K*dt),'*','linewidth',4)
hold on
plot(T,alpha*S(:,2),'--','linewidth',2)
 
hold off
 xlabel('\alpha t','fontsize',20)
ylabel('dp_0/dt','fontsize',20)
set(gca, 'XTick', 0:2:14);
axis([0 14 0 0.4])
 

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

 

 
 
    