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

# include "fsolve.h"
# include "b1g3.h"

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

void b1g3 ( void dydt ( double t, double y[], double dy[] ), 
  double tspan[2], double y0[], int n, int m, double t[], double y[] )

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

    b1g3() uses the b1g3 method + fsolve_be() to solve an ODE.

  Licensing:

    This code is distributed under the MIT license.

  Modified:

    28 November 2025

  Author:

    John Burkardt

  Input:

    dydt: a function that evaluates the right hand side of the ODE,
    of the form
      void dydt ( double t, double y[], double dy[] )

    double tspan[2]: the initial and final times.

    double y0[m]: 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 dt;
  double *fm;
  int i;
  int info;
  int j;
  double t1;
  double t2;
  double t3;
  double tol;
  double *y1;
  double *y2;
  double *y3;
  double *yp1;
  double *yp2;

  fm = ( double * ) malloc ( m * sizeof ( double ) );
  y1 = ( double * ) malloc ( m * sizeof ( double ) );
  y2 = ( double * ) malloc ( m * sizeof ( double ) );
  y3 = ( double * ) malloc ( m * sizeof ( double ) );
  yp1 = ( double * ) malloc ( m * sizeof ( double ) );
  yp2 = ( double * ) malloc ( m * sizeof ( double ) );

  dt = ( tspan[1] - tspan[0] ) / ( double ) ( n );

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

  info = 0;
  tol = 1.0e-05;

  for ( i = 0; i <= n; i++ )
  {
/*
  Step 0 uses the initial condition.
*/
    if ( i == 0 )
    {
      for ( j = 0; j < m; j++ )
      {
        y[i*m+j] = y0[j];
      }
    }
/*
  Step 1 uses the midpoint method.
*/
    else if ( i == 1 )
    {
      t1 = t[i];
      for ( j = 0; j < m; j++ )
      {
        y1[j] = y[(i-1)*m+j];
      }
      dydt ( t1, y1, yp1 );

      t2 = t1 + 0.5 * dt;
      for ( j = 0; j < m; j++ )
      {
        y2[j] = y1[j] + 0.5 * dt * yp1[j];
      }

      info = fsolve_be ( dydt, m, t1, y1, t2, y2, fm, tol );

      if ( info != 1 )
      {
        printf ( "\n" );
        printf ( "b1g3(): Fatal error in fsolve_be()!\n" );
        printf ( "  info = %d\n", info );
        exit ( 1 );
      }

      for ( j = 0; j < m; j++ )
      {
        y[i*m+j] = 2.0 * y2[j] - y1[j];
      }
    }
/*
  Later steps use B1G3.
*/  
    else
    {
      t1 = t[i-2];
      for ( j = 0; j < m; j++ )
      {
        y1[j] = y[(i-2)*m+j];
      }

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

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

      info = fsolve_b1g3 ( dydt, m, t1, y1, t2, y2, t3, y3, fm, tol );

      if ( info != 1 )
      {
        printf ( "\n" );
        printf ( "b1g3(): Fatal error in fsolve_be()!\n" );
        printf ( "  info = %d\n", info );
        exit ( 1 );
      }

      for ( j = 0; j < m; j++ )
      {
        y[i*m+j] = y3[j];
      }
    }
  }
/*
  Free memory.
*/
  free ( fm );
  free ( y1 );
  free ( y2 );
  free ( y3 );
  free ( yp1 );
  free ( yp2 );

  return;
}