/* Runge-Kutta-4 integrator x(t), scalar version. */

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

float dxdt(float x, float t)
{
    return x;
}

float true_solution(float x0, float t)

/* Modify this when dxdt is changed! */

{
    return x0 * exp(t);
}

void output(float x0, float x, float t)
{
    printf("%f %f %f\n", t, x, x - true_solution(x0, t));
}

main()
{
    /* Fourth-order Runge-Kutta integrator x(t). */

    float x, x0, t, dt;

    /* Initial conditions: */

    t = 0;
    x = x0 = 1;

    dt = .1;

    while (t < 1) {

	float dx1, dx2, dx3, dx4;

	output(x0, x, t);

	dx1 = dt*dxdt(x, t);
	dx2 = dt*dxdt(x + 0.5*dx1, t + 0.5*dt);
	dx3 = dt*dxdt(x + 0.5*dx2, t + 0.5*dt);
	dx4 = dt*dxdt(x + dx3, t + dt);

	t += dt;
	x += (dx1 + 2*dx2 + 2*dx3 + dx4)/6.0;
    }
    output(x0, x, t);
}