Function finite_diff_jacobian 

pub fn finite_diff_jacobian<S: Scalar>(
    r: &dyn Fn(&[S], &mut [S]),
    x: &[S],
    m: usize,
    jac: &mut [S],
)

Compute a finite-difference Jacobian of residual r at x using the canonical central FD step cbrt(S::EPSILON) * (1 + |x|), which scales with |x| to remain valid for large-magnitude parameters. The additive scaling prevents silent zero-Jacobian failures that occur when an unscaled step falls below the precision floor of x + h == x (manifesting for |x| > ~5e7 at f64 precision).

jac is row-major m*n.