/*
* euler.c
* (c) Neil Gershenfeld  3/2/03
* integrate simple harmonic motion with an Euler algorithm
*/

#include <math.h>

double error,final_err,t,h,tmax,pi,y[2],ynew[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) {
		ynew[0] = y[0] + h * y[1];
		ynew[1] = y[1] - h * y[0];
		y[0] = ynew[0];
		y[1] = ynew[1];
		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]);
}