function [ S, I, R ] = sir ( m, beta, gamma )

%% sir follows an SIR model for m steps.
%
%  Input, integer M, the number of steps to take.  1000 might be typical.
%
%  Input, real BETA, GAMMA, the infection and recovery rates.
%  0.02 and 0.005 might be typical.
%
%  Output, real S(M), I(M), R(M), the susceptible, infected, and recovered
%  populations at each time step.
%

%
%  Allocate space.
%
  S = zeros ( 1, m );
  I = zeros ( 1, m );
  R = zeros ( 1, m );

  for i = 1 : m

    if ( i == 1 )
      S(i) = 99;
      I(i) = 1;
      R(i) = 0;
      N = S(i) + I(i) + R(i);
    else
      S(i) = S(i-1) - beta * S(i-1) * I(i-1) / N;
      I(i) = I(i-1) + beta * S(i-1) * I(i-1) / N - gamma * I(i-1);
      R(i) = R(i-1)                              + gamma * I(i-1);
    end

  end
  
  return
end
