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

  global N Y0

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

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

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