% Gillespie algorithm for a simple 2 state model

% initialize the states
clear
N = 100;

alpha = 2;
beta = 1;


K = 100; % number of steps
T = zeros(N,1);
P = zeros(N,1);

for j = 2:K
    % pick a random number
    indx = mod(j,2);
    k = rand(N,1);
    stepper = (1-indx)*alpha+indx*beta;
    T(:,j) = -log(k)/stepper+T(:,j-1);
    P(:,j) = mod(j,2)*ones(N,1);
   
end

plot(T',P')

%prepare stsatistics
%find the smallest value of t
klast = K;
t_min = min(T(:,klast));
dt = t_min/K;

for k = 1:K-1
 for j = 1:N
     tindx = min(find(T(j,:)>=k*dt));
     p(j) = P(j,tindx);
 end
 pt(k) = sum(p)/N;
tm(k) = k*dt;

sol(k) = alpha/(alpha+beta)*(1-exp(-(alpha+beta)*tm(k)));
end
figure(2)
sig = sqrt((sol-sol.^2)/N)
plot(tm,pt,tm,sol,tm,sol+sig,'--',tm,sol-sig,'--')

axis([0 50 0 1])