
%%%%% SIR stochastic simulation

sizes=zeros(1000,1);
for j=1:1000   % do 1000 realizations of the model

t=0;
n=10000;  %% code is slow for n=10000 !!
i=1;   % one initial infective

s=n-i;
tc=i;

beta=100.0;
gamma=50.0;

beta_over_n=beta/n;

while (i>0)
    r(1)=beta_over_n*s*i; % infection rate
    r(2)=gamma*i; % recovery rate
    
    total_rate=r(1)+r(2);
        
     
    t=t+exprnd(1.0/total_rate);
    
    x=total_rate*rand(1);
    if (x<r(1))    
           s=s-1;
           i=i+1;
           tc=tc+1;    
       else       
           i=i-1;
                  
    end
    
    
end
j,tc  % print out realization number and outcome
      % (mainly so that we can see our progress)
sizes(j)=tc;
end
hist(sizes,100)  
% the "100" means use 100 bars in the histogram


    