/*
* rk.c
* (c) Neil Gershenfeld  3/2/03
* integrate simple harmonic motion
*    with 4th order Runge-Kutta
*/

#include <math.h>

double error,final_err,t,h,tmax,pi,y[2],ynew[2],
   A[2],B[2],C[2],D[2];
int nsteps;

main() {
	printf("Step size?  "); scanf("%lf",&h);
	pi = 4.0 * atan(1.0);
	tmax = 100*pi;
	t = 0;
	nsteps = ((int) (tmax / h));
	y[0] = 1;
	y[1] = 0;
	error = 0;
	while (t < tmax) {
		A[0] = h * y[1];
		A[1] = -h * y[0];
		B[0] = h * (y[1] + A[1]/2.0);
		B[1] = -h * (y[0] + A[0]/2.0);
		C[0] = h * (y[1] + B[1]/2.0);
		C[1] = -h * (y[0] + B[0]/2.0);
		D[0] = h * (y[1] + C[1]);
		D[1] = -h * (y[0] + C[0]);
		y[0] = y[0] + (A[0]+2.0*B[0]+2.0*C[0]+D[0])/6.0;
		y[1] = y[1] + (A[1]+2.0*B[1]+2.0*C[1]+D[1])/6.0;
		t += h;
		error = error + fabs(cos(t) - y[0]);
	}
	error = error / nsteps;
	final_err = fabs(cos(t) - y[0]);
	printf("error: %g  final error: %g  slope: %g\n",error, final_err, y[1]);
}