diff --git a/src/main.rs b/src/main.rs
index 191fac4..fd30a2c 100644
--- a/src/main.rs
+++ b/src/main.rs
@@ -3,6 +3,7 @@ extern crate nalgebra as na;
 
 mod univariate_solvers;
 mod univariate_minimizers;
+mod simulated_annealing;
 mod nelder_mead;
 mod particle_swarm_optimization;
 mod non_linear_least_squares;
@@ -131,7 +132,7 @@ fn test_univariate_solvers(verbose: bool) {
     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();
-    let x_itp :       f64 = univariate_solvers::itp_solve(&(fct as fn(f64) -> f64),  -5.0, 1.0, 0.1, 2.0, 1.0, tol);
+    let x_itp :       f64 = univariate_solvers::itp_solve(&(fct as fn(f64) -> f64),  -5.0, 1.0, 0.1, 2.0, 1.0, tol).unwrap();
     num_tests_passed += check_result(x_newton, x_mathematica, tol, "Newton's method", verbose); num_tests_total += 1;
     num_tests_passed += check_result(x_newton_num, x_mathematica, tol, "Newton's method (num)", verbose); num_tests_total += 1;
     num_tests_passed += check_result(x_halley, x_mathematica_2, tol, "Halley's method", verbose); num_tests_total += 1;
@@ -145,6 +146,7 @@ fn test_univariate_solvers(verbose: bool) {
 }
 
 fn test_univariate_optimizers(verbose: bool) {
+    println!("Testing univariate optimizers.");
     let tol :      f64 = 1e-10;
     // let max_iter : u32 = 100;
     // let dx_num :   f64 = 1e-6;
@@ -159,6 +161,7 @@ fn test_univariate_optimizers(verbose: bool) {
 }
 
 fn test_multivariate_optimizers(verbose: bool) {
+    println!("Testing multivariate optimizers.");
     let tol :      f64 = 1e-10;
     let max_iter : u32 = 1000;
     // let dx_num :   f64 = 1e-6;
@@ -168,24 +171,26 @@ fn test_multivariate_optimizers(verbose: bool) {
     let tol_x:      f64 = 1e-4;
     let tol_f_x:    f64 = 1e-6;
     
-    let x_true:        na::DVector<f64> = na::DVector::from_vec(vec![1.,1.]);
-    let f_x_true:      f64 = rosenbrock(&x_true);
-    let sol_nelder_mead: (na::DVector<f64>, f64) = nelder_mead::nelder_mead_minimize(&(rosenbrock as fn(&na::DVector<f64>) -> f64), &na::DVector::from_vec(vec![2.0,-1.0]), 0.1, tol, max_iter, false);
-    // num_tests_passed += check_result(x_nelder_mead, x_true, tol*1e2, "Nelder-Mead", verbose);
+    let x_init: na::DVector<f64> = na::DVector::from_vec(vec![2.0,-1.0]);
+    let x_true: na::DVector<f64> = na::DVector::from_vec(vec![1.,1.]);
+    let f_x_true: f64 = rosenbrock(&x_true);
+    let sol_nelder_mead: (na::DVector<f64>, f64) = nelder_mead::nelder_mead_minimize(&(rosenbrock as fn(&na::DVector<f64>) -> f64), &x_init, 0.1, tol, max_iter, false);
+    let sol_simulated_annealing: (na::DVector<f64>, f64) = simulated_annealing::simulated_annealing(&rosenbrock, &x_init, 0.1, 100.0, 0.95, max_iter, 123);
     num_tests_passed += check_result_optim(&sol_nelder_mead.0, sol_nelder_mead.1, &x_true, f_x_true, tol_x, tol_f_x, "Nelder-Mead", verbose); num_tests_total += 1;
+    num_tests_passed += check_result_optim(&sol_simulated_annealing.0, sol_simulated_annealing.1, &x_true, f_x_true, 0.1, 0.01, "Simulated Annealing", verbose); num_tests_total += 1;
     print_test_results(num_tests_passed, num_tests_total);
 }
 
-fn test_particle_swarm_debug() {
-    let n_particles: u32 = 10;
-    let lb: na::DVector<f64> = na::DVector::from_vec(vec![-5.0,-5.0]);
-    let ub: na::DVector<f64> = na::DVector::from_vec(vec![5.0,5.0]);
-    let tol: f64 = 1e-6;
-    let n_iter_max: u32 = 1000;
-    let rng_seed: u32 = 69;
-    let x: f64 = particle_swarm_optimization::particle_swarm_minimize(&rosenbrock, n_particles, &lb, &ub, tol, n_iter_max, rng_seed);
-    println!("x = {}", x);
-}
+// fn test_particle_swarm_debug() {
+//     let n_particles: u32 = 10;
+//     let lb: na::DVector<f64> = na::DVector::from_vec(vec![-5.0,-5.0]);
+//     let ub: na::DVector<f64> = na::DVector::from_vec(vec![5.0,5.0]);
+//     let tol: f64 = 1e-6;
+//     let n_iter_max: u32 = 1000;
+//     let rng_seed: u32 = 69;
+//     let x: f64 = particle_swarm_optimization::particle_swarm_minimize(&rosenbrock, n_particles, &lb, &ub, tol, n_iter_max, rng_seed);
+//     println!("x = {}", x);
+// }
 
 fn test_non_linear_lsqr_solvers(verbose: bool) {
     let mut num_tests_passed : u32 = 0;
diff --git a/src/simulated_annealing.rs b/src/simulated_annealing.rs
new file mode 100644
index 0000000..fb4a692
--- /dev/null
+++ b/src/simulated_annealing.rs
@@ -0,0 +1,95 @@
+extern crate nalgebra as na;
+
+#[path = "./xorwow.rs"]
+mod xorwow;
+
+use xorwow::Xorwow;
+use na::DVector;
+
+/// Performs simulated annealing optimization on a given objective function.
+///
+/// Simulated annealing is a probabilistic technique for approximating the global optimum of a given function.
+/// It is inspired by the annealing process in metallurgy, where a material is heated and then slowly cooled
+/// to increase the size of its crystals and reduce defects.
+///
+/// # Arguments
+///
+/// * `f` - A closure that takes a vector of `f64` values and returns a single `f64` value representing the objective function to be optimized.
+/// * `initial_x` - The initial guess or starting point for the optimization process, represented as a `DVector<f64>`.
+/// * `search_scale` - A scaling factor that determines the magnitude of the random perturbations applied during the optimization process.
+/// * `temperature_0` - The initial temperature for the simulated annealing process.
+/// * `temperature_decay` - The rate at which the temperature decreases during the optimization process.
+/// * `num_iterations` - The maximum number of iterations to perform during the optimization process.
+///
+/// # Returns
+///
+/// The function returns a `DVector<f64>` representing the approximate global optimum of the objective function.
+///
+/// # Example
+///
+/// ```
+/// extern crate nalgebra as na;
+/// 
+/// let rosenbrock = |x: &na::DVector<f64>| -> f64 {
+///     return (1.0-x[0]).powi(2) + 100.0*(x[1] - x[0].powi(2)).powi(2);
+/// };
+/// 
+/// // Initial guess
+/// let initial_x: na::DVector<f64> = na::DVector::from_vec(vec![0.0, 0.0]);
+/// 
+/// // Simulated annealing parameters
+/// let search_scale = 0.1;
+/// let temperature_0 = 100.0;
+/// let temperature_decay = 0.99;
+/// let num_iterations = 1000;
+/// let seed = 123;
+/// 
+/// // Perform simulated annealing optimization
+/// let result = simulated_annealing::simulated_annealing(&rosenbrock, &initial_x, search_scale, temperature_0, temperature_decay, num_iterations, seed);
+/// 
+/// println!("Approximate global optimum: {:?}", result);
+/// ```
+pub fn simulated_annealing(
+    f: &dyn Fn(&DVector<f64>) -> f64,
+    initial_x: &DVector<f64>,
+    search_scale: f64,
+    temperature_0: f64,
+    temperature_decay: f64,
+    num_iterations: u32,
+    seed: u32
+) -> (na::DVector<f64>, f64) {
+    let mut temperature = temperature_0;
+    let mut rng = Xorwow::new(seed);
+    let mut current_x = initial_x.clone();
+    let mut best_x = current_x.clone();
+    let mut current_energy = f(&current_x);
+    let mut best_energy = current_energy;
+    
+    for _ in 0..num_iterations {
+        // Generate a random vector
+        let mut next_x = current_x.clone();
+        for i in 0..current_x.len() {
+            next_x[i] += search_scale * (2.0 * rng.next_f64() - 1.0);
+        }
+
+        // Calculate the energy of the new vector
+        let next_energy = f(&next_x);
+        
+        if next_energy < best_energy {
+            best_energy = next_energy;
+            best_x = next_x.clone();
+        }
+
+        // Decide whether to accept the new vector based on Metropolis-Hastings criterion
+        if next_energy < current_energy
+            || f64::exp((current_energy - next_energy) / temperature) > rng.next_f64()
+        {
+            current_x = next_x.clone();
+            current_energy = next_energy;
+        }
+
+        temperature *= temperature_decay;
+    }
+    
+    (best_x, best_energy)
+}
diff --git a/src/univariate_solvers.rs b/src/univariate_solvers.rs
index a9049f8..8ad0cee 100644
--- a/src/univariate_solvers.rs
+++ b/src/univariate_solvers.rs
@@ -233,6 +233,15 @@ fn swap(a: &mut f64, b: &mut f64) {
     *b = temp;
 }
 
+/// @brief Solves the equation f(x) = 0 for x in [a, b]. The root must be bracketed by [a, b], meaning that f(a) * f(b) < 0.
+/// @param `f` A reference to a function that takes a `f64` as an argument and returns a `f64`.
+/// @param `a` The left endpoint of the interval [a, b].
+/// @param `b` The right endpoint of the interval [a, b].
+/// @param `k1` Hyperparameter in [0, inf]
+/// @param `k2` Hyperparameter in [1, 1 + phi], with phi = (1 + sqrt(5))/2
+/// @param `n0` Hyperparameter in [0, inf]
+/// @param `epsilon` Tolerance on x.
+/// @return Solution if f(a) * f(b) < 0, an error string otherwise.
 pub fn itp_solve(
     f: &dyn Fn(f64) -> f64,
     a: f64,
@@ -241,7 +250,7 @@ pub fn itp_solve(
     k2: f64,
     n0: f64,
     epsilon: f64
-) -> f64 {
+) -> Result<f64, &'static str> {
     let mut a: f64 = a;
     let mut b: f64 = b;
 
@@ -253,6 +262,10 @@ pub fn itp_solve(
         swap(&mut a, &mut b);
         swap(&mut ya, &mut yb);
     }
+
+    if ya * yb > 0.0 {
+        return Err("Root is not bracketed")
+    }
     
     // Preprocessing
     let nhalf: f64 = f64::log2((b - a)/(2.0 * epsilon));
@@ -301,5 +314,5 @@ pub fn itp_solve(
         j += 1;
     }
 
-    (a + b) / 2.0
+    Ok((a + b) / 2.0)
 }