% Gillespie algorithm for a simple 2 state model

% initialize the states
clear

%set the rate constants
c(1) = 1;
c(2) = 0.005;
c(3) = 0.6;

%Specify the change matrix

C = [1,0;-1,1;0,-1];

%initialize the state space

x(1) = 15;
x(2) = 10;

X = x;
T(1) = 0;
j = 1;
H = 1;
K = 10000; % number of time steps
while (H>0&j<K)
    j = j+1;
    h(1) = c(1)*x(1);
    h(2) = c(2)*x(1)*x(2);
    h(3) = c(3)*x(2);
    hc = cumsum(h);
    H = sum(h);
    % if H = 0, the population went extinct
    rn = rand(1);
    if (H>0)
    T(j) = -log(rn)/H +T(j-1);

    rn = rand(1);
    rk = min(find(rn<=hc/H)); % this determines which reaction occurs
    x = x + C(rk,:);
    X(j,:) = x;
    end
end

figure(1)
plot(T,X)

figure(2)
plot(X(:,1),X(:,2))