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

  global N

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

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