let y' t y = cos t
let exact t = sin t

let rk4_step y t h =
  let k1 = h *. y' t y in
  let k2 = h *. y' (t +. 0.5*.h) (y +. 0.5*.k1) in
  let k3 = h *. y' (t +. 0.5*.h) (y +. 0.5*.k2) in
  let k4 = h *. y' (t +. h) (y +. k3) in
  y +. (k1+.k4)/.6.0 +. (k2+.k3)/.3.0

let rec loop steps h n y t =
  if n < steps then loop steps h (n+1) (rk4_step y t h) (t +. h) else
    y
let _ =
  let h = 0.1 in
  let y = loop 102 h 1 1.0 0.0 in
  let err = abs_float (y -. (exact ((float_of_int 102) *. h))) in
  let large = 10000000 in
  let y = loop large h 1 1.0 0.0 in
  Printf.printf "%b\n"
    (abs_float (y -. (exact (float_of_int large) *. h)) < 2. *. err)