#include "iostream.h"
#include "math.h"

ab4am4(double h, double xstart, double y, double xend)
{
	double x = xstart, k1, k2, k3, k4;
	double yy[4], yp;
	yy[0] = y;
	for(int i = 1; i < 4; i++)
	{
		k1 = feuler(x, yy[i - 1]);
		k2 = feuler(x + h / 2, yy[i - 1] + h * k1 / 2);
		k3 = feuler(x + h / 2, yy[i - 1] + h * k2 / 2);
		k4 = feuler(x + h, yy[i - 1] + h * k3);
		yy[i] = yy[i - 1] + h * (k1 + 2 * k2 + 2 * k3 + k4) / 6;

		x += h;
		double ytemp = 1 / (1 + x * x);
		cout << i << "\t" << x << "\t" << yy[i] << "\t" << ytemp << "\t" << sqrt((yy[i] - ytemp)*(yy[i] - ytemp)) << endl;
	}
	for(int j = 4; j < (xend-xstart) / h; j++)
	{
		yp = yy[3] + h / 24 * (55 * feuler(x, yy[3]) - 59 * feuler(x - h, yy[2])
			+ 37 * feuler(x - 2 * h, yy[1]) - 9 * feuler(x - 3 * h, yy[0]));
		y = yy[3] + h / 24 * (9 * feuler(x + h, yp) + 19 * feuler(x, yy[3])
			- 5 * feuler(x - h, yy[2]) + feuler(x - 2 * h, yy[1]));
//		y = (y + yp) / 2;

		x += h;
		yy[0] = yy[1];
		yy[1] = yy[2];
		yy[2] = yy[3];
		yy[3] = y;
		double ytemp = 1 / (1 + x * x);
		cout << j << "\t" << x << "\t" << y << "\t" << ytemp << "\t" << sqrt((y - ytemp)*(y - ytemp)) << endl;
	}

	return 0;	
}