% Order: S, E, ES, P
clear AllTs AllNs
N = [1000;1000;0;0];
Tmax = 2;
index=1;
V=100;
AllNs(index,:)=N;
AllTs(index)=0;
rf=1/V;
rr=1;
rcat=1;
Results = [-1 1 0; -1 1 1; 1 -1 -1; 0 0 1];
time=0;
while time < Tmax
    % Calculate rates
    Rates = [rf*N(2)*N(1) rr*N(3) rcat*N(3)];
    NextTimes = -log(rand(1,3))./Rates;
    [tmin,kmin]=min(NextTimes);
    time = time+tmin;
    if (time < Tmax)
        N=N+Results(:,kmin);
        index=index+1;
        AllTs(index)=time;
        AllNs(index,:)=N;
    end
end