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

  global N 

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

%

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