function sir_continuous ( t0, tmax, k_user, tau_user )

  global k
  global tau

  k = k_user;
  tau = tau_user;

  tspan = [ t0, tmax ];

  y0 = [ 0.995, 0.005, 0.000 ];

  [ tout, yout ] = ode23 ( @epidemic, tspan, y0 );

  plot ( tout, yout(:,1), 'g', ...
         tout, yout(:,2), 'r', ...
         tout, yout(:,3), 'k' );

  return
end
function dydt = epidemic ( t, y )

  global k
  global tau

% tau = 0.8;
% k = 4;

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

  dydt(1,1) = - tau * I * S;
  dydt(2,1) =   tau * I * S - I / k;
  dydt(3,1) =                 I / k;

  return
end
