function stochastic_SIR_closed_five_realizations
% store and plot five realizations of the model
% note that our plot respects the "jump process" nature of this model, so
% we have to store two values of prevalence at each time point (the
% previous and the new value)

number_of_realizations=5;
prev_results=NaN(2000,number_of_realizations);  % array where we will store prevalences
time_results=NaN(2000,number_of_realizations);  % array where we will store corresponding times
% we use NaN so that unwanted points don't get plotted
for loop=1:number_of_realizations
    
t=0; N=1000;
I=1; S=N-I;  % one initial infective, everyone else susceptible
beta=100.0; gamma=50.0;

prev_results(1,loop)=I;
time_results(1,loop)=t;
output_index=2;


while (I > 0)      % simulation will end when infection goes extinct

    r(1)=beta*S*I/N; % infection rate
    r(2)=gamma*I; % recovery rate   

    total_rate=r(1)+r(2);

    t=t+exprnd(1.0/total_rate);  % generates exponential random variable

    prev_results(output_index,loop)=I;
    time_results(output_index,loop)=t;
% looks odd, but because we have a jump process, I remains at its previous
% value until time t, then jumps up or down by one

    x=total_rate*rand(1);   % generates uniform r.v. on (0,total_rate)

    if (x<r(1))    
           S=S-1;
           I=I+1;
	else       
           I=I-1;
    end
    
    prev_results(output_index+1,loop)=I;
    time_results(output_index+1,loop)=t;
    output_index=output_index+2;
    
end

    plot(time_results(:,1),prev_results(:,1),...
        time_results(:,2),prev_results(:,2),...
        time_results(:,3),prev_results(:,3),...
        time_results(:,4),prev_results(:,4),...
        time_results(:,5),prev_results(:,5) )
    

    
end

