# include <cstdlib>
# include <iostream>
# include <cmath>

using namespace std;

double midpoint ( double t0, double u0, double dt, 
  double f ( double t, double u ) );
double f3 ( double t, double u );

int main ( )
//
//  MIDPOINT_F3 uses the midpoint method (more accurate than Euler)
//  on problem F3, the "wiggly" function.
//
{
  double dt, exact, pi = 3.14159265, t0, t1, tmax, u0, u1;

  u0 = 0.5;
  t0 = 0.0;
  tmax = 12.0 * pi;
  
  dt = 0.25;
    
  while ( true )
  {
    exact = 0.5 * exp ( sin ( t0 ) );

    cout << "  " << t0 
         << "  " << u0 
         << "  " << exact 
         << "  " << fabs ( u0 - exact ) << "\n";

    if ( tmax <= t0 )
    {
      break;
    }
    t1 = t0 + dt;
    u1 = midpoint ( t0, u0, dt, f3 );

    t0 = t1;
    u0 = u1;
  }
  
  return 0;
}
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;
}  

double f3 ( double t, double u )
{
  double dudt;
  
  dudt = u * cos ( t );
  
  return dudt;
}