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

using namespace std;

void f99 ( double t, double u[], double dudt[] );

int main ( )
//
//  MIDPOINT_F99 uses the midpoint method on problem F7, the pendulum function.
//
{
  double dt, dudt[2], t0, t1, th, tmax, u0[2], u1[2], uh[2];
  int i;

  u0[0] = 0.0;
  u0[1] = 0.0;
  t0 = 0.0;
  tmax = 20.0;
  
  dt = 0.04;
    
  while ( true )
  {
    cout << "  " << t0 
         << "  " << u0[0] 
         << "  " << u0[1] << "\n";

    if ( tmax <= t0 )
    {
      break;
    }
//
//  Half step.
//
    f99 ( t0, u0, dudt );
    th = t0 + 0.5 * dt;
    for ( i = 0; i < 2; i++ )
    {
      uh[i] = u0[i] + 0.5 * dt * dudt[i];
    }
//
//  Full step.
//
    f99 ( th, uh, dudt );
    t1 = t0 + dt;
    for ( i = 0; i < 2; i++ )
    {
      u1[i] = u0[i] + dt * dudt[i];
    }
//
//  Shift.
//
    t0 = t1;
    for ( i = 0; i < 2; i++ )
    {
      u0[i] = u1[i];
    }
  }
  
  return 0;
}
void f99 ( double t, double u[], double dudt[] )
{
  dudt[0] = u[1];
  dudt[1] = pow ( u[0], 3 ) / 6.0 - u[0] + 2.0 * sin ( 2.7853 * t );

  return;
}