Function minimize_bfgs 

pub fn minimize_bfgs(
    f: &dyn Fn(&[f64]) -> f64,
    grad: &dyn Fn(&[f64]) -> Vec<f64>,
    x0: &[f64],
    tol: f64,
    max_iter: usize,
) -> OptResult

BFGS quasi-Newton method with Armijo line search.

Maintains an approximate inverse Hessian H (initialized to identity) and updates it with the BFGS rank-2 formula at each step.