//   Runge Kutta algorithm for first-order differential equations

//   To adjust for new cases 
//   you have to redefine the function f(x,y) as in 
//	y' = f(x,y)
//   see double f () at end of file.

//   The program has to request from the user also :

//      Initial value for x and y , use variables x0 and y0
//      the upper bound for x
//	the step size

//
//
//  Compile + link with    g++ runge_kutta.cpp -o runge_kutta
//
//  Program uses  math-functions , user-defined functions,
//
//


#include <iostream>
#include <cmath>
#include <iomanip>    // Manipulating output formats
using namespace std;

double runge4( double , double , double );	// Runge-Kutta function
double f     ( double , double ) ;              // as y'=f(x,y)

/// Insert your main function at this spot
//

///////////   end of main  /////////////////////////////
//
// 4th order Runge-Kutta method , error scales with stepsize to the 4th power

double runge4 (double x, double y, double step) {

   double h, k1, k2, k3, k4;

   h = step / 2.0 ;
   k1=step*f(x, y);
   k2=step*f(x+h, y+k1/2.0);
   k3=step*f(x+h, y+k2/2.0);
   k4=step*f(x+step, y+k3);

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

/////////////////////////////////////////////////
//  Test function :
//  y' = - y  solution is A*e^{-x} , starting with x=0 and y=1 we have A=1.0
//  with e^{-4} = 0.018315639

double  f (double x, double y)
{
  return(-y);
}