/// @brief Newton's method for solving a function f(x) = 0
/// @param f function to solve
/// @param df derivative of function f
/// @param x0 initial guess
/// @param tol tolerance
/// @param max_iter maximum number of iterations
/// @return solution
pub fn newton_solve<F: Fn(f64) -> f64, F2: Fn(f64) -> f64>(f : &F, df : &F2, x0 : f64, tol : f64, max_iter : u32) -> f64 {
    let mut x: f64 = x0;
    let mut dx: f64;
    let mut fx: f64;
    let mut dfx: f64;
    for _iter in 0..max_iter {
        fx = f(x);
        dfx = df(x);
        if dfx == 0.0 {
            dx = fx;
        } else {
            dx = fx/dfx;
        }
        x -= dx;
        if f64::abs(dx) < tol {
            break;
        }
    }
    return x;
}

/// @brief Newton's method for solving a function f(x) = 0
/// @param f function to solve
/// @param df derivative of function f
/// @param x0 initial guess
/// @param tol tolerance
/// @param max_iter maximum number of iterations
/// @return solution
pub fn newton_solve_num<F: Fn(f64) -> f64>(f : &F, x0 : f64, tol : f64, dx_num : f64, max_iter : u32) -> f64 {
    return newton_solve(f, &((|x: f64| (f(x + dx_num) - f(x - dx_num))/(2.0*dx_num)) as fn(f64) -> f64), x0, tol, max_iter);
}