clear all
clear reset
 
% This program is written by kevin Burrage (UQ) and Margherita Carletti (Urbino) 
% define parameters  
% The rate constants are in the vector c

A=6.0221415*10^23;
V=10^(-15);
AV=A*V;
c=[0.078 0.043 0.0007 0.0039 0.05*10^9/AV 0.5 0.012*10^9/AV 0.4791 0.00012*10^9/AV 0.8765*10^(-11)]
t_end = 100;
X_initial=[2 6 0 2 0 0]'

% use the ssa
[X,t] = ssa_goutsias(c,t_end,X_initial);

plot(t, X(1,1:end), t, X(2,1:end), t, X(3,1:end), t, X(4,1:end),t, X(5,1:end),t, X(6,1:end))
xlabel('time')
ylabel('number')
legend('m','D','RNA','DNA','DNA1','DNA2',1)
title('SSA solver')