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

  global N

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

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