%
%              hiv_rhc
%

  function dy = hiv_rhs(t,y,params)

  global coef Y0 wt 

  bE = params(1);
  delta = params(2);
  d1 = params(3);
  k2 = params(4);
  lam1 = params(5);
  Kb = params(6);

  epsilon = coef(1);
  k1 = coef(2);
  lam2 = coef(3);
  d2 = coef(4);
  f = coef(5); 
  NT = coef(6);
  c = coef(7);
  rho1 = coef(8);
  rho2 = coef(9);
  lamE = coef(10);
  dE = coef(11);  
  deltaE = coef(12); 
  m1 = coef(13);
  m2 = coef(14);
  Kd = coef(15);

  dy = [lam1 - d1*y(1) - (1 - epsilon)*k1*y(5)*y(1);
        lam2 - d2*y(2) - (1 - f*epsilon)*k2*y(5)*y(2);
        (1 - epsilon)*k1*y(5)*y(1) - delta*y(3) - m1*y(6)*y(3);
        (1 - f*epsilon)*k2*y(5)*y(2) - delta*y(4) - m2*y(6)*y(4);
        NT*delta*(y(3)+y(4)) - c*y(5) - ((1 - epsilon)*rho1*k1*y(1) + (1 - f*epsilon)*rho2*k2*y(2))*y(5);
        lamE + (bE*(y(3)+y(4))./(y(3)+y(4)+Kb))*y(6) - (dE*(y(3)+y(4))./(y(3)+y(4)+Kd))*y(6) - deltaE*y(6)];