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

  global N 

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

  S = y(1);
  I = y(2);
  S_delta = y(3);
  I_delta = y(4);
  S_gamma = y(5);
  I_gamma = y(6);
  S_r = y(7);
  I_r = y(8);
  Svec = [S_delta; I_delta; S_gamma; I_gamma; S_r; I_r];
  J = [(-delta-gamma*I) -gamma*S; gamma*I (gamma*S-(r+delta))];
  Der = blkdiag(J,J,J); 
  Grad = [(N-S); -I; -I*S; I*S; 0; -I];
  Sen = Der*Svec + Grad;

%
% Construct the right-hand side
%
  dy = [delta*(N-S) - gamma*I*S;
        gamma*I*S - (r+delta)*I;
        Sen];