%
%          SIR_rhs
%
  function dy = SIR_rhs(t,y,params);

  global N Y0

  gamma = params(1);
  r = params(2);
  delta = params(3);
  k = params(4);

%
% Exploit the fact that R(t) = N - S(t) - I(t).
%

  dy = [delta*N - delta*y(1) - gamma*k*y(2)*y(1);
        gamma*k*y(2)*y(1) - (r + delta)*y(2)];