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

  global N

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

  S = y(1);
  I = y(2);
  R = y(3);

  dy = [mu*(N-S) - eta*k*I*S;
        eta*k*I*S - (gamma + mu)*I;
        gamma*I - mu*R];