The Runge-Kutta method (in its midpoint incarnation) for a system of differential equations can be implemented, e.g. like

K0 = F(x, Y0); K1/2 = F(x + h/2; Y0 + h/2 K0); Y1 = Y0 + hK1/2; err = max(|K0 - K1/2|) h/2

The Predictor-Corrector method can be generalized as, e.g.

Prediction:

Y(x) ≈ Y0+F0⋅(x-x0)+C⋅(x-x0)2, where C=[Y-1-Y0-F0⋅(x-1-x0)]/(x-1-x1)2 <= from Y(x-1)=Y-1;
Ypred=Y(x1); F1=F(x1,Ypred)

Correction:

Y(x) ≈ Y0+F0⋅(x-x0)+C⋅(x-x0)2 +D⋅(x-x0)2(x-x-1),
where D=(F1-F0-2C⋅(x1-x0))/(2(x1-x0)(x1-x-1)+(x1-x0)2) <= from Y'(x1)=F1;
Ycorr=Y(x1)

Adaptive step-size algorithm, written in (supposedly selfexplanatory) Haskell, recursive:

runge_kutta:: RealFloat z =>
          (z->[z]->z) -> [z] -> [[z]] -> z -> z -> z -> z -> ([z],[z])

runge_kutta   f           x       y      b    h    eps  acc
	| last x >= b 			= (x,y)
	| (last x)+h > b 		= runge_kutta f x     y     b last_h eps acc
	| err <= 2*tol			= runge_kutta f new_x new_y b new_h  eps acc
	| otherwise			= runge_kutta f x     y     b new_h  eps acc
	where
		last_h = b - last x
		tol = (eps*abs(y1)+acc)*sqrt(abs( h/(b - a) ))
		a = head x
		new_x = x ++ [last x + h]
		new_y = y ++ [y1]
		(y1,err) = stepper_routine f (last x) (last y) h
		new_h = h*(tol/err)**0.2*0.9

Problems

1. Generalize your routines to system of differential equations. Let your driver estimate the global_error=√(∑_i err_i²).
2. Integrate the equation y''(x)=-y , y(0)=0 , y'(0)=1 from a=0 to b=(4+¹/₂)π.