Function central_difference 

pub fn central_difference<F>(f: F, x: f64, h: f64) -> f64
where
    F: Fn(f64) -> f64,

Approximates the first derivative with a central difference.

This computes (f(x + h) - f(x - h)) / (2h). When h == 0.0, the result may be infinite or NaN according to normal floating-point behavior.