/* calculation of integral rk4 */
#include <iostream>
#include <iomanip>
#include <cmath>
using namespace std;


int main()
{
    double t = 0.;
    double y = 1.;
    double h = 0.001;
    double k1,k2,k3,k4;

    while(t<1.) {
        k1  = y;
        k2  = y+0.5*h*k1;
        k3  = y+0.5*h*k2;
        k4  = y+h*k3;
        y  += h/6.*(k1+2.*k2+2.*k3+k4);
        t  += h;
    }

    double y_exact = exp(t);

    cout << "RK4: " << setprecision(18) << y << " Exact: " << setprecision(18) << y_exact <<  \
    " diff: " << fabs(y-y_exact) << endl;


//    printf("Euler: %.16f, exact: %.16f, diff: %.16f\n",y,y_exact,fabs(y-y_exact));
}