double midpoint ( double t0, double u0, double dt, 
  double f ( double t, double u ) )
{
  double th, uh, u1;
//
//  Take a half time step and estimate UH there.
//
  th = t0 + 0.5 * dt;
  uh = u0 + 0.5 * dt * f ( t0, u0 );
//
//  Evaluate the derivative at (TH,UH).
//
  u1 = u0 +       dt * f ( th, uh );
  
  return u1;
}