diff --git a/src/main.rs b/src/main.rs
index f8974c9..817196e 100644
--- a/src/main.rs
+++ b/src/main.rs
@@ -39,6 +39,7 @@ fn main() {
     let x_newton = univariate_solvers::newton_solve(&(fct as fn(f64) -> f64), &(dfct as fn(f64) -> f64), x0, tol, max_iter);
     let x_newton_num: f64 = univariate_solvers::newton_solve_num(&(fct as fn(f64) -> f64), x0, tol, dx_num, max_iter);
     let x_halley: f64 = univariate_solvers::halley_solve(&(fct as fn(f64) -> f64), &(dfct as fn(f64) -> f64), &(ddfct as fn(f64) -> f64), x0, tol, max_iter, false).unwrap();
+    let x_halley_num: f64 = univariate_solvers::halley_solve_num(&(fct as fn(f64) -> f64), x0, tol, dx_num, max_iter, false).unwrap();
     let x_bisection : f64 = univariate_solvers::bisection_solve(&(fct as fn(f64) -> f64), -5.0, 1.0, tol).unwrap();
     let x_secant : f64 = univariate_solvers::secant_solve(&(fct as fn(f64) -> f64), -1.0, 1.0, tol, max_iter);
     let x_ridder : f64 = univariate_solvers::ridder_solve(&(fct as fn(f64) -> f64), -5.0, 1.0, tol, max_iter).unwrap();
@@ -46,6 +47,7 @@ fn main() {
     println!("Newton's method       : x = {}\tf(x) = {}", x_newton, fct(x_newton));
     println!("Newton's method (num) : x = {}\tf(x) = {}", x_newton_num, fct(x_newton_num));
     println!("Halley's method       : x = {}\tf(x) = {}", x_halley, fct(x_halley));
+    println!("Halley's method (num) : x = {}\tf(x) = {}", x_halley_num, fct(x_halley_num));
     println!("Bisection             : x = {}\tf(x) = {}", x_bisection, fct(x_bisection));
     println!("Secant method         : x = {}\tf(x) = {}", x_secant, fct(x_secant));
     println!("Ridder's method       : x = {}\tf(x) = {}", x_ridder, fct(x_ridder));
diff --git a/src/univariate_solvers.rs b/src/univariate_solvers.rs
index 69bc98d..4a100ea 100644
--- a/src/univariate_solvers.rs
+++ b/src/univariate_solvers.rs
@@ -30,11 +30,12 @@ pub fn newton_solve<F, F2>(f : F, df : F2, x0 : f64, tol : f64, max_iter : u32)
 
 /// @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 dx_num numerical differentiation step size
 /// @param max_iter maximum number of iterations
 /// @return solution
+/// @note This method uses numerical differentiation to compute the first and second derivatives.
 pub fn newton_solve_num<F>(f : F, x0 : f64, tol : f64, dx_num : f64, max_iter : u32) -> f64
 where F : Fn(f64) -> f64
 {
@@ -74,6 +75,41 @@ where F : Fn(f64) -> f64, F2 : Fn(f64) -> f64, F3 : Fn(f64) -> f64
     return Err("Halley method did not converge after reaching the maximum number of iterations allowed.")
 }
 
+/// Halley's method for solving a function f(x) = 0
+/// @param f function to solve
+/// @param x0 initial guess
+/// @param tol tolerance
+/// @param dx_num numerical differentiation step size
+/// @param max_iter maximum number of iterations
+/// @return solution
+/// @note This method is more efficient than Newton's method, but requires the second derivative of f.
+/// @note This method uses numerical differentiation to compute the first and second derivatives.
+pub fn halley_solve_num<F>(f: F, x0: f64, tol: f64, dx_num : f64, max_iter: u32, verbose: bool) -> Result<f64, &'static str>
+where F : Fn(f64) -> f64
+{
+    let mut x: f64 = x0;
+    let mut f_x: f64;
+    let mut df_x: f64;
+    let mut ddf_x: f64;
+    let mut fx_m_dx: f64;
+    let mut fx_p_dx: f64;
+    for _i in 0..max_iter {
+        f_x = f(x);
+        fx_m_dx = f(x - dx_num);
+        fx_p_dx = f(x + dx_num);
+        df_x = (fx_p_dx - fx_m_dx)/(2.0*dx_num);
+        ddf_x = (fx_p_dx - 2.0*f_x + fx_m_dx) / (dx_num.powi(2));
+        if verbose {
+            println!("x = {}, f(x) = {}, df(x) = {}, ddf(x) = {}", x, f_x, df_x, ddf_x);
+        }
+        if f64::abs(f_x) < tol {
+            return Ok(x);
+        }
+        x = x - 2.0*f_x*df_x / (2.0*df_x.powi(2) - f_x*ddf_x);
+    }
+    return Err("Halley method did not converge after reaching the maximum number of iterations allowed.")
+}
+
 // --------------------------------------------------------------------
 // ------------------------ Bracketing methods ------------------------
 // --------------------------------------------------------------------