# include <stdbool.h>
# include <stdio.h>
# include <stdlib.h>

# include "bdf3.h"
# include "fsolve.h"

/******************************************************************************/

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

/******************************************************************************/
/*
  Purpose:

    bdf3() approximates the solution to an ODE using the bdf3 method.

  Licensing:

    This code is distributed under the MIT license.

  Modified:

    29 May 2025

  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 *k1;
  double *k2;
  double *k3;
  double t1;
  double t2;
  double t3;
  double t4;
  double tol;
  double *y1;
  double *y2;
  double *y3;
  double *y4;

  dely = ( double * ) malloc ( m * sizeof ( double ) );
  fn = ( double * ) malloc ( m * sizeof ( double ) );
  k1 = ( double * ) malloc ( m * sizeof ( double ) );
  k2 = ( double * ) malloc ( m * sizeof ( double ) );
  k3 = ( double * ) malloc ( m * sizeof ( double ) );
  y1 = ( double * ) malloc ( m * sizeof ( double ) );
  y2 = ( double * ) malloc ( m * sizeof ( double ) );
  y3 = ( double * ) malloc ( m * sizeof ( double ) );
  y4 = ( double * ) malloc ( m * sizeof ( double ) );

  for ( i = 0; i < n + 1; i++ )
  {
    t[i] = ( ( double ) ( n - i ) * tspan[0] 
           + ( double ) (     i ) * tspan[1] ) 
           / ( double ) ( n     );
  }

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

  for ( i = 0; i <= n; i++ )
  {
/*
  Initial condition.
*/
    if ( i == 0 )
    {
      for ( j = 0; j < m; j++ )
      {
        y[i*m+j] = y0[j];
      }
    }
/*
  Explicit Runge-Kutta order 3
*/
    else if ( i < 3 )
    {
      for ( j = 0; j < m; j++ )
      {
        y1[j] = y[(i-1)*m+j];
      }
      dydt ( t[i-1], y1, k1 );

      for ( j = 0; j < m; j++ )
      {
        y2[j] = y[(i-1)*m+j] + dt * k1[j];
      }
      dydt ( t[i-1] + dt, y2, k2 );

      for ( j = 0; j < m; j++ )
      {
        y3[j] = y[(i-1)*m+j] + 0.25 * dt * ( k1[j] + k2[j] );
      }
      dydt ( t[i-1] + 0.5 * dt, y3, k3 );

      for ( j = 0; j < m; j++ )
      {
        y[i*m+j] = y[(i-1)*m+j] + dt * ( k1[j] + k2[j] + 4.0 * k3[j] ) / 6.0;
      }
    }
/*
  BDF3
*/
    else
    {
      t1 = t[i-3];
      t2 = t[i-2];
      t3 = t[i-1];
      t4 = t[i];
      dydt ( t[i-1], y+(i-3)*m, dely );
      for ( j = 0; j < m; j++ )
      {
        y4[j] = y[(i-1)*m+j] + dt * dely[j];
      }

      info = fsolve_bdf3 ( dydt, m, t1, y+(i-3)*m, t2, y+(i-2)*m, t3, 
        y+(i-1)*m, t4, y4, fn, tol );

      if ( info != 1 )
      {
        printf ( "\n" );
        printf ( "bdf3(): Fatal error!\n" );

        printf ( "  info = %d\n", info );

        exit ( 1 );
      }

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

  free ( dely );
  free ( fn );
  free ( k1 );
  free ( k2 );
  free ( k3 );
  free ( y1 );
  free ( y2 );
  free ( y3 );
  free ( y4 );

  return;
}