let simpson f a b n =
  let h = (b -. a) /. float_of_int n in
  let sum = ref (f a +. f b) in
  for i = 1 to n - 1 do
    let x = a +. float_of_int i *. h in
    let coef = if i mod 2 = 0 then 2.0 else 4.0 in
    sum := !sum +. coef *. f x
  done;
  h *. !sum /. 3.0

;;

int_of_float (simpson (fun x -> x *. x) 0.0 1.0 100 *. 30000.0)