# include <cstdlib>
# include <iomanip>
# include <iostream>

using namespace std;

# include "bdf2.hpp"
# include "fsolve.hpp"

//****************************************************************************80

void bdf2 ( void dydt ( double x, double y[], double yp[] ), double tspan[2], 
  double y0[], int n, int m, double t[], double y[] )

//****************************************************************************80
//
//  Purpose:
//
//    bdf2() approximates the solution to an ODE using the BDF2 method.
//
//  Licensing:
//
//    This code is distributed under the MIT license.
//
//  Modified:
//
//    01 December 2023
//
//  Author:
//
//    John Burkardt
//
//  Input:
//
//    dydt: a function that evaluates the right hand side of the ODE.
//
//    double tspan[2]: contains the initial and final times.
//
//    double y0[m]: a column vector containing the initial condition.
//
//    int n: the number of steps to take.
//
//    int m: the number of variables.
//
//  Output:
//
//    double t[n+1], y[m*(n+1)]: the times and solution values.
//
{
  double *dely;
  double dt;
  double *fn;
  int i;
  int info;
  int j;
  double t1;
  double t2;
  double t3;
  double tol;
  double *y1;
  double *y2;
  double *y3;

  dely = new double[m];
  fn = new double[m];
  y1 = new double[m];
  y2 = new double[m];
  y3 = new double[m];

  dt = ( tspan[1] - tspan[0] ) / ( double ) ( n );
  tol = 1.0E-05;

  for ( i = 0; i <= n; i++ )
  {
    if ( i == 0 )
    {
      t[i] = tspan[0];
      for ( j = 0; j < m; j++ )
      {
        y[i*m+j] = y0[j];
      }
    }
    else if ( i == 1 )
    {
      t1 = t[i-1];
      for ( j = 0; j < m; j++ )
      {
        y1[j] = y[(i-1)*m+j];
      }

      t2 = t1 + dt;
      dydt ( t1, y1, dely );
      for ( j = 0; j < m; j++ )
      {
        y2[j] = y1[j] + dt * dely[j];
      }
      info = fsolve_be ( dydt, m, t1, y1, t2, y2, fn, tol );
      if ( info != 1 )
      {
        cout << "\n";
        cout << "bdf2(): Fatal error!\n";
        cout << "  info = " << info << "\n";
        exit ( 1 );
      }

      t[i] = t2;
      for ( j = 0; j < m; j++ )
      {
        y[i*m+j] = y2[j];
      }
    }
    else
    {
      t1 = t[i-2];
      for ( j = 0; j < m; j++ )
      {
        y1[j] = y[(i-2)*m+j];
      }
      t2 = t1 + dt;
      for ( j = 0; j < m; j++ )
      {
        y2[j] = y[(i-1)*m+j];
      }
      t3 = t2 + dt;

      dydt ( t2, y2, dely );
      for ( j = 0; j < m; j++ )
      {
        y3[j] = y2[j] + dt * dely[j];
      }

      info = fsolve_bdf2 ( dydt, m, t1, y1, t2, y2, t3, y3, fn, tol );

      if ( info != 1 )
      {
        cout << "\n";
        cout << "bdf2(): Fatal error!\n";
        cout << "  info = " << info << "\n";
        exit ( 1 );
      }

      t[i] = t3;
      for ( j = 0; j < m; j++ )
      {
        y[i*m+j] = y3[j];
      }
    }
  }

  delete [] dely;
  delete [] fn;
  delete [] y1;
  delete [] y2;
  delete [] y3;

  return;
}