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

  global N 

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

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

%
%
  dy = [delta*(N-S) - gamma*I*S;
        gamma*I*S - (r+delta)*I];