Numerical integration using adaptive Simpson's rule with substitution of variables to handle infinite limits

NumericalIntegration.tq

let integrate (x0, x1) f =
  let integrateAux f (x0, x2) =
    let δ = ε in
    let simpsonsRule dx fx0 fx1 fx2 = dx * (fx0 + 4*fx1 + fx2) / 6 in
    let rec aux x0 fx0 x2 fx2 x4 fx4 =
      let dx = x2 - x0 in
      let x1 = (x0 + x2) / 2 in
      let x3 = (x2 + x4) / 2 in
      if x0 = x1 || x2 = x3 then 0 else
        let fx1 = f x1 in
        let fx3 = f x3 in
        let y = simpsonsRule (2 * dx) fx0 fx2 fx4 in
        let v = max (abs x0) (abs x4) in
        let i1 = simpsonsRule dx fx0 fx1 fx2 in
        let i2 = simpsonsRule dx fx2 fx3 fx4 in
        if dx < δ * v || abs(i1 + i2 - y) < δ then y else
          aux x0 fx0 x1 fx1 x2 fx2 + aux x2 fx2 x3 fx3 x4 fx4 in
    let x1 = (x2 + x0)/2 in
    aux x0 (f x0) x1 (f x1) x2 (f x2) in

  let handleInfiniteLimits x0 x1 f =
    let x0, x1, f = if x0 ≤ x1 then x0, x1, f else x1, x0, [x → -f x] in
    let δ = ε in
    let g1 t = f(x1 - (1 - t) / t) / (t*t) in
    let g2 t = f(x0 + t / (1 - t)) / sqr(1 - t) in
    let g3 t = (1 + t*t) / sqr(1 - t*t) * f(t / (1 - t*t)) in
    (x0 = -∞, x1 = ∞) @
    [ False, False → x0, x1, f
    | True, False → δ, 1-δ, g1
    | False, True → δ, 1-δ, g2
    | True, True → δ-1, 1-δ, g3 ] in
  
  let x0, x1, f = handleInfiniteLimits x0 x1 f in
  
  integrateAux f (x0, x1)