/* Runge-Kutta method of order 4 for first order ODEs 
 *
 * Implemented in ANSI C.
 *
 * Author: Michael Mazack, <mazack @ yahoo . com>
 *
 * Date: March 29th, 2006
 *
 * License: Public Domain. Redistribution and modification without 
 * restriction is granted. 
 */
#include <stdio.h>

/* Prototype for our ODE, dx/dt = f(t, x) */
double ode1(double, double);

/* Prototype for Runge-Kutta method of order 4 */
double rk4(double, double, double, double (*)(double, double));

int main()
{
  /* Initial conditions:
   * x(0) = 1 
   */
  double t = 0;
  double x = 1;

  /* Step size */
  double h = 0.001;

  /* Here we find all RK4 estimates on the interval [0,1] for t in
   * intervals of h. 
   */
  for(t=0;t<=1;t+=h)
    {
      x = rk4(t, x, h, &ode1);
      printf("%lf\t%lf\n", (t + h), x);
    }

  return 0;
}

/* Note: Any first order ODE can be used here.
 * This function is only a sample.
 */
double ode1(double t, double x)
{
  /* dx/dt = -3*x^2 + 2*t */
  return (-3*x*x + 2*t);
}

/* Implementation of the Runge-Kutta method of order 4 for
 * estimating differential equations.
 *
 * First argument: t (independent variable of the differential equation)
 * Second argument: x (initial condition or previous evaluation of rk4())
 * Third argument: h (step size)
 * Fourth argument: (*f) (a function pointer to dx/dt = f(t, x)) 
 */
double rk4(double t, double x, double h, double (*f)(double, double))
{
  /* Since this function will be called many many times, it
   * is probably best to declare these as static. 
   */
  static double k1, k2, k3, k4;

  /* This is one of the most elegant things I've ever seen. */ 
  k1 = h*f(t, x);
  k2 = h*f(t + 0.5*h, x + 0.5*k1);
  k3 = h*f(t + 0.5*h, x + 0.5*k2);
  k4 = h*f(t + h, x + k3);

  return (x + (1.0/6.0)*(k1 + 2*k2 + 2*k3 + k4));
}