// Midpoint method y(x), modular form, array version.

#include <iostream>
#include <cmath>
using namespace std;

#define real double
#define N 2

void get_dydx(real y[], real x, real dydx[])	// specify the ODE
{
    dydx[0] = y[1];
    dydx[1] = -y[0];
}

void step(real y[], real& x, real dx)		// Midpoint integrator
{
    int i;
    real dydx[N], dy[N], y1[N];

    get_dydx(y, x, dydx);

    for (i = 0; i < N; i++) {
	dy[i] = dx*dydx[i];
	y1[i] = y[i] + 0.5*dy[i];
    }

    get_dydx(y1, x+0.5*dx, dydx);

    for (i = 0; i < N; i++) {
	dy[i] = dx*dydx[i];
	y[i] += dy[i];
    }

    x += dx;
}

bool time_to_stop(real y[], real x, real x_max)
{
    return (x >= x_max);
}

void output(real y[], real x)
{
    int i;

    printf("%f ", x);
    for (i = 0; i < N; i++)
	printf("%f ", y[i]);
    printf("\n");
}

void initialize(real y[], real& x, real& dx, real& x_max)
{
    // Specify initial conditions.

    x = 0;

    y[0] = 0;
    y[1] = 1;

    dx = .01;
    x_max = 1;
}

main()
{
    real y[N], x, dx, x_max;

    initialize(y, x, dx, x_max);

    while (!time_to_stop(y, x, x_max)) {
	output(y, x);
	step(y, x, dx);
    }

    output(y, x);
}