# include <cmath>
# include <cstdlib>

using namespace std;

# include "rk3.hpp"

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

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

//****************************************************************************80
//
//  Purpose:
//
//    rk3() uses a Runge-Kutta order 3 explicit method to solve an ODE.
//
//  Licensing:
//
//    This code is distributed under the MIT license.
//
//  Modified:
//
//    07 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 dt;
  int i;
  int j;
  double *k1;
  double *k2;
  double *k3;
  double *ym;

  ym = new double[m];

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

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

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

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

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

    t[j] = t[j-1] + dt;
    for ( i = 0; i < m; i++ )
    {
      y[i+j*m] = y[i+(j-1)*m] + dt * ( k1[i] + k2[i] + 4.0 * k3[i] ) / 6.0;
    }
    delete [] k1;
    delete [] k2;
    delete [] k3;
  }

  delete [] ym;

  return;
}