double rk4( double x0, double y0, double h, double (*fun)(double, double))
{
	double y1;
	double k1,k2,k3,k4;
	k1 = h*fun(x0     ,y0      );
	k2 = h*fun(x0+h/2.,y0+k1/2.);
	k3 = h*fun(x0+h/2.,y0+k2/2.);
	k4 = h*fun(x0+h   ,y0+k3   );
	y1 = y0 + (k1+2.*k2+2.*k3+k4)/6.;
	return y1;
}