% Linda J. S. Allen TTU
%CTMC SIR Epidemic with No Births and Deaths
%Three Sample paths and Time Epidemic Ended
clear all
beta=0.5;  g=.25; N=100;
time=50; init=2; R0=beta/g
for k=1:3 
    clear  t s i
    t(1)=0;  i(1)=init; s(1)=N-i(1);
    ss(1)=s(1); ii(1)=i(1);
    j=1;
    %Gillespie's Algorithm
    while i(j)>0 & j<=1000
        y1=rand; y2=rand;
        t(j+1)=-log(y1)/((beta/N)*i(j)*s(j)+g*i(j))+t(j);
        if (y2<=(beta/N)*s(j)/(beta/N*s(j)+g)) 
	       i(j+1)=i(j)+1;
           s(j+1)=s(j)-1;
        else
           i(j+1)=i(j)-1;
           s(j+1)=s(j);   
       end
       j=j+1;
    end
    time_end(k)=t(j); % Time  epidemic ends.
    if k==1
       l1=stairs(t,i,'r'); set(l1,'LineWidth',2);
       hold on
    elseif k==2
        l1=stairs(t,i,'g'); set(l1,'LineWidth',2);
    else
        l1=stairs(t,i,'b'); set(l1,'LineWidth',2);
    end
end
time_end %Time epidemic ends
dt=.01; %Euler's method for ODEs
for tt=1:time/dt
   ss(tt+1)=ss(tt)+dt*(-beta*ss(tt)*ii(tt)/N);
   ii(tt+1)=ii(tt)+dt*(beta*ss(tt)*ii(tt)/N -g*ii(tt));
end
set(gca,'fontsize',18);
plot([0:dt:time],ii,'k--','LineWidth',2); 
ylabel('I(t)'); xlabel('t'); title('CTMC SIR Model')
hold off