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

using namespace std;

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

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

  u0[0] = pi / 3;
  u0[1] = 0.0;
  t0 = 0.0;
  tmax = pi;
  
  dt = 0.025;
    
  while ( true )
  {
    cout << "  " << t0 
         << "  " << u0[0] 
         << "  " << u0[1] << "\n";

    if ( tmax <= t0 )
    {
      break;
    }
//
//  Half step.
//
    f7 ( 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.
//
    f7 ( 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 f7 ( double t, double u[], double dudt[] )
{
  double g = 386.09, l = 10.0;

  dudt[0] = u[1];
  dudt[1] = - ( g / l ) * sin ( u[0] );

  return;
}