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

  global N 

  beta = params(1);
  gamma = params(2);
  mu = params(3);

  S = y(1);
  I = y(2);
  S_S0 = y(3);
  I_S0 = y(4);
  S_I0 = y(5);
  I_I0 = y(6);
  Svec = [S_S0; I_S0; S_I0; I_I0];
  J = [(-mu-beta*I) -beta*S; beta*I (beta*S-(gamma+mu))];
  Der = blkdiag(J,J); 
  Sen = Der*Svec;

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