diff --git a/Cargo.lock b/Cargo.lock
index d3386e9..8abff31 100644
--- a/Cargo.lock
+++ b/Cargo.lock
@@ -187,11 +187,21 @@ version = "1.13.0"
 source = "registry+https://github.com/rust-lang/crates.io-index"
 checksum = "60b1af1c220855b6ceac025d3f6ecdd2b7c4894bfe9cd9bda4fbb4bc7c0d4cf0"
 
+[[package]]
+name = "float-cmp"
+version = "0.9.0"
+source = "registry+https://github.com/rust-lang/crates.io-index"
+checksum = "98de4bbd547a563b716d8dfa9aad1cb19bfab00f4fa09a6a4ed21dbcf44ce9c4"
+dependencies = [
+ "num-traits",
+]
+
 [[package]]
 name = "ganesh"
 version = "0.9.0"
 dependencies = [
  "criterion",
+ "float-cmp",
  "nalgebra",
  "num",
 ]
diff --git a/Cargo.toml b/Cargo.toml
index f59ff3e..a8feb63 100644
--- a/Cargo.toml
+++ b/Cargo.toml
@@ -18,6 +18,7 @@ num = "0.4.3"
 
 [dev-dependencies]
 criterion = "0.5.1"
+float-cmp = "0.9.0"
 
 [[bench]]
 name = "rosenbrock_benchmark"
diff --git a/README.md b/README.md
index fabe10a..617392f 100644
--- a/README.md
+++ b/README.md
@@ -65,21 +65,31 @@ To minimize this function, we could consider using the Nelder-Mead algorithm:
 use ganesh::prelude::*;
 use ganesh::algorithms::NelderMead;
 
-let problem = Rosenbrock { n: 2 };
-let nm = NelderMead::default();
-let mut m = Minimizer::new(nm, 2);
-let x0 = &[2.0, 2.0];
-m.minimize(&problem, x0, &mut ()).unwrap();
-println!("{}", m.status);
+fn main() -> Result<(), Infallible> {
+    let problem = Rosenbrock { n: 2 };
+    let nm = NelderMead::default();
+    let mut m = Minimizer::new(nm, 2);
+    let x0 = &[2.0, 2.0];
+    m.minimize(&problem, x0, &mut ())?;
+    println!("{}", m.status);
+    Ok(())
+}
 ```
 
 This should output
 ```shell
 MSG:       term_f = STDDEV
-X:         [0.9999999946231828, 0.9999999884539057]
-F(X):      0.00000000000000009170942877687133
-N_EVALS:   160
+X:         +1.000 ± 0.707
+           +1.000 ± 1.416
+F(X):      +0.000
+N_F_EVALS: 159
+N_G_EVALS: 0
 CONVERGED: true
+COV:       
+  ┌             ┐
+  │ 0.500 1.000 │
+  │ 1.000 2.005 │
+  └             ┘
 ```
 
 # Bounds
diff --git a/src/algorithms/bfgs.rs b/src/algorithms/bfgs.rs
index 5d9e376..4d373fb 100644
--- a/src/algorithms/bfgs.rs
+++ b/src/algorithms/bfgs.rs
@@ -1,7 +1,7 @@
 use nalgebra::{DMatrix, DVector, RealField, Scalar};
 use num::Float;
 
-use crate::{Algorithm, Bound, Function, Status};
+use crate::{convert, Algorithm, Bound, Function, Status};
 
 use super::line_search::{LineSearch, StrongWolfeLineSearch};
 
@@ -68,12 +68,12 @@ pub struct BFGS<T: Scalar, U, E> {
     status: Status<T>,
     x: DVector<T>,
     g: DVector<T>,
-    p: DVector<T>,
     h_inv: DMatrix<T>,
     f_previous: T,
     terminator_f: BFGSFTerminator<T>,
     terminator_g: BFGSGTerminator<T>,
     line_search: Box<dyn LineSearch<T, U, E>>,
+    max_step: T,
     error_mode: BFGSErrorMode,
 }
 
@@ -117,7 +117,6 @@ where
             status: Default::default(),
             x: Default::default(),
             g: Default::default(),
-            p: Default::default(),
             h_inv: Default::default(),
             f_previous: T::infinity(),
             terminator_f: BFGSFTerminator {
@@ -127,6 +126,7 @@ where
                 tol_g_abs: Float::cbrt(T::epsilon()),
             },
             line_search: Box::new(StrongWolfeLineSearch::default()),
+            max_step: convert!(1e8, T),
             error_mode: Default::default(),
         }
     }
@@ -150,7 +150,7 @@ where
 
 impl<T, U, E> Algorithm<T, U, E> for BFGS<T, U, E>
 where
-    T: RealField + Float,
+    T: RealField + Float + Default,
 {
     fn initialize(
         &mut self,
@@ -159,6 +159,8 @@ where
         bounds: Option<&Vec<Bound<T>>>,
         user_data: &mut U,
     ) -> Result<(), E> {
+        self.status = Status::default();
+        self.f_previous = T::infinity();
         self.h_inv = DMatrix::identity(x0.len(), x0.len());
         self.x = Bound::to_unbounded(x0, bounds);
         self.g = func.gradient_bounded(self.x.as_slice(), bounds, user_data)?;
@@ -178,23 +180,23 @@ where
         bounds: Option<&Vec<Bound<T>>>,
         user_data: &mut U,
     ) -> Result<(), E> {
-        self.p = -&self.h_inv * &self.g;
+        let d = -&self.h_inv * &self.g;
         let (valid, alpha, f_kp1, g_kp1) = self.line_search.search(
             &self.x,
-            &self.p,
-            None,
+            &d,
+            Some(self.max_step),
             func,
             bounds,
             user_data,
             &mut self.status,
         )?;
         if valid {
-            let s = self.p.scale(alpha);
+            let dx = d.scale(alpha);
             let grad_kp1_vec = g_kp1;
-            let y = &grad_kp1_vec - &self.g;
+            let dg = &grad_kp1_vec - &self.g;
             let n = self.x.len();
-            self.update_h_inv(i_step, n, &s, &y);
-            self.x += s;
+            self.update_h_inv(i_step, n, &dx, &dg);
+            self.x += dx;
             self.g = grad_kp1_vec;
             self.status
                 .update_position((Bound::to_bounded(self.x.as_slice(), bounds), f_kp1));
@@ -236,12 +238,49 @@ where
                 if covariance.is_none() {
                     covariance = hessian.pseudo_inverse(Float::cbrt(T::epsilon())).ok();
                 }
-                self.status.cov = covariance
+                self.status.set_cov(covariance);
             }
             BFGSErrorMode::ApproximateHessian => {
-                self.status.cov = Some(self.h_inv.clone());
+                self.status.set_cov(Some(self.h_inv.clone()));
             }
         }
         Ok(())
     }
 }
+
+#[cfg(test)]
+mod tests {
+    use std::convert::Infallible;
+
+    use float_cmp::approx_eq;
+
+    use crate::{prelude::*, test_functions::Rosenbrock};
+
+    use super::BFGS;
+
+    #[test]
+    fn test_bfgs() -> Result<(), Infallible> {
+        let algo = BFGS::default();
+        let mut m = Minimizer::new(algo, 2).with_max_steps(10000);
+        let problem = Rosenbrock { n: 2 };
+        m.minimize(&problem, &[-2.0, 2.0], &mut ())?;
+        assert!(m.status.converged);
+        assert!(approx_eq!(f64, m.status.fx, 0.0, epsilon = 1e-10));
+        m.minimize(&problem, &[2.0, 2.0], &mut ())?;
+        assert!(m.status.converged);
+        assert!(approx_eq!(f64, m.status.fx, 0.0, epsilon = 1e-10));
+        m.minimize(&problem, &[2.0, -2.0], &mut ())?;
+        assert!(m.status.converged);
+        assert!(approx_eq!(f64, m.status.fx, 0.0, epsilon = 1e-10));
+        m.minimize(&problem, &[-2.0, -2.0], &mut ())?;
+        assert!(m.status.converged);
+        assert!(approx_eq!(f64, m.status.fx, 0.0, epsilon = 1e-10));
+        m.minimize(&problem, &[0.0, 0.0], &mut ())?;
+        assert!(m.status.converged);
+        assert!(approx_eq!(f64, m.status.fx, 0.0, epsilon = 1e-10));
+        m.minimize(&problem, &[1.0, 1.0], &mut ())?;
+        assert!(m.status.converged);
+        assert!(approx_eq!(f64, m.status.fx, 0.0, epsilon = 1e-10));
+        Ok(())
+    }
+}
diff --git a/src/algorithms/lbfgs.rs b/src/algorithms/lbfgs.rs
index b53f853..213db7b 100644
--- a/src/algorithms/lbfgs.rs
+++ b/src/algorithms/lbfgs.rs
@@ -3,7 +3,7 @@ use std::collections::VecDeque;
 use nalgebra::{DVector, RealField, Scalar};
 use num::Float;
 
-use crate::{Algorithm, Bound, Function, Status};
+use crate::{convert, Algorithm, Bound, Function, Status};
 
 use super::line_search::{LineSearch, StrongWolfeLineSearch};
 
@@ -68,7 +68,6 @@ pub struct LBFGS<T: Scalar, U, E> {
     status: Status<T>,
     x: DVector<T>,
     g: DVector<T>,
-    p: DVector<T>,
     f_previous: T,
     terminator_f: LBFGSFTerminator<T>,
     terminator_g: LBFGSGTerminator<T>,
@@ -76,6 +75,7 @@ pub struct LBFGS<T: Scalar, U, E> {
     m: usize,
     y_store: VecDeque<DVector<T>>,
     s_store: VecDeque<DVector<T>>,
+    max_step: T,
     error_mode: LBFGSErrorMode,
 }
 
@@ -124,7 +124,6 @@ where
             status: Default::default(),
             x: Default::default(),
             g: Default::default(),
-            p: Default::default(),
             f_previous: T::infinity(),
             terminator_f: LBFGSFTerminator {
                 tol_f_abs: T::epsilon(),
@@ -136,6 +135,7 @@ where
             m: 10,
             y_store: VecDeque::default(),
             s_store: VecDeque::default(),
+            max_step: convert!(1e8, T),
             error_mode: Default::default(),
         }
     }
@@ -148,7 +148,7 @@ where
     fn g_approx(&self) -> DVector<T> {
         let m = self.s_store.len();
         let mut q = self.g.clone();
-        if m <= 2 {
+        if m < 1 {
             return q;
         }
         let mut a = vec![T::zero(); m];
@@ -171,7 +171,7 @@ where
 
 impl<T, U, E> Algorithm<T, U, E> for LBFGS<T, U, E>
 where
-    T: RealField + Float,
+    T: RealField + Float + Default,
 {
     fn initialize(
         &mut self,
@@ -180,6 +180,8 @@ where
         bounds: Option<&Vec<Bound<T>>>,
         user_data: &mut U,
     ) -> Result<(), E> {
+        self.status = Status::default();
+        self.f_previous = T::infinity();
         self.x = Bound::to_unbounded(x0, bounds);
         self.g = func.gradient_bounded(self.x.as_slice(), bounds, user_data)?;
         self.status.inc_n_g_evals();
@@ -198,34 +200,38 @@ where
         bounds: Option<&Vec<Bound<T>>>,
         user_data: &mut U,
     ) -> Result<(), E> {
-        self.p = -self.g_approx();
+        let d = -self.g_approx();
         let (valid, alpha, f_kp1, g_kp1) = self.line_search.search(
             &self.x,
-            &self.p,
-            None,
+            &d,
+            Some(self.max_step),
             func,
             bounds,
             user_data,
             &mut self.status,
         )?;
         if valid {
-            let dx = self.p.scale(alpha);
-            self.s_store.push_back(dx.clone());
-            if self.s_store.len() > self.m {
-                self.s_store.pop_front();
-            }
+            let dx = d.scale(alpha);
             let grad_kp1_vec = g_kp1;
-            self.y_store.push_back(&grad_kp1_vec - &self.g);
-            if self.y_store.len() > self.m {
-                self.y_store.pop_front();
+            let dg = &grad_kp1_vec - &self.g;
+            let sy = dx.dot(&dg);
+            let yy = dg.dot(&dg);
+            if sy > T::epsilon() * yy {
+                self.s_store.push_back(dx.clone());
+                self.y_store.push_back(dg);
+                if self.s_store.len() > self.m {
+                    self.s_store.pop_front();
+                    self.y_store.pop_front();
+                }
             }
             self.x += dx;
             self.g = grad_kp1_vec;
             self.status
                 .update_position((Bound::to_bounded(self.x.as_slice(), bounds), f_kp1));
         } else {
-            self.status.update_message("LINE SEARCH FAILED");
-            self.status.set_converged();
+            // reboot
+            self.s_store.clear();
+            self.y_store.clear();
         }
         Ok(())
     }
@@ -262,9 +268,46 @@ where
                 if covariance.is_none() {
                     covariance = hessian.pseudo_inverse(Float::cbrt(T::epsilon())).ok();
                 }
-                self.status.cov = covariance
+                self.status.set_cov(covariance);
             }
         }
         Ok(())
     }
 }
+
+#[cfg(test)]
+mod tests {
+    use std::convert::Infallible;
+
+    use float_cmp::approx_eq;
+
+    use crate::{prelude::*, test_functions::Rosenbrock};
+
+    use super::LBFGS;
+
+    #[test]
+    fn test_lbfgs() -> Result<(), Infallible> {
+        let algo = LBFGS::default();
+        let mut m = Minimizer::new(algo, 2);
+        let problem = Rosenbrock { n: 2 };
+        m.minimize(&problem, &[-2.0, 2.0], &mut ())?;
+        assert!(m.status.converged);
+        assert!(approx_eq!(f64, m.status.fx, 0.0, epsilon = 1e-10));
+        m.minimize(&problem, &[2.0, 2.0], &mut ())?;
+        assert!(m.status.converged);
+        assert!(approx_eq!(f64, m.status.fx, 0.0, epsilon = 1e-10));
+        m.minimize(&problem, &[2.0, -2.0], &mut ())?;
+        assert!(m.status.converged);
+        assert!(approx_eq!(f64, m.status.fx, 0.0, epsilon = 1e-10));
+        m.minimize(&problem, &[-2.0, -2.0], &mut ())?;
+        assert!(m.status.converged);
+        assert!(approx_eq!(f64, m.status.fx, 0.0, epsilon = 1e-10));
+        m.minimize(&problem, &[0.0, 0.0], &mut ())?;
+        assert!(m.status.converged);
+        assert!(approx_eq!(f64, m.status.fx, 0.0, epsilon = 1e-10));
+        m.minimize(&problem, &[1.0, 1.0], &mut ())?;
+        assert!(m.status.converged);
+        assert!(approx_eq!(f64, m.status.fx, 0.0, epsilon = 1e-10));
+        Ok(())
+    }
+}
diff --git a/src/algorithms/lbfgsb.rs b/src/algorithms/lbfgsb.rs
index e28fd63..4586884 100644
--- a/src/algorithms/lbfgsb.rs
+++ b/src/algorithms/lbfgsb.rs
@@ -368,7 +368,7 @@ where
 
 impl<T, U, E> Algorithm<T, U, E> for LBFGSB<T, U, E>
 where
-    T: RealField + Float + TotalOrder,
+    T: RealField + Float + TotalOrder + Default,
 {
     fn initialize(
         &mut self,
@@ -377,6 +377,9 @@ where
         bounds: Option<&Vec<Bound<T>>>,
         user_data: &mut U,
     ) -> Result<(), E> {
+        self.status = Status::default();
+        self.f_previous = T::infinity();
+        self.theta = T::one();
         self.l = DVector::from_element(x0.len(), T::neg_infinity());
         self.u = DVector::from_element(x0.len(), T::infinity());
         if let Some(bounds_vec) = bounds {
@@ -496,9 +499,46 @@ where
                 if covariance.is_none() {
                     covariance = hessian.pseudo_inverse(Float::cbrt(T::epsilon())).ok();
                 }
-                self.status.cov = covariance
+                self.status.set_cov(covariance);
             }
         }
         Ok(())
     }
 }
+
+#[cfg(test)]
+mod tests {
+    use std::convert::Infallible;
+
+    use float_cmp::approx_eq;
+
+    use crate::{prelude::*, test_functions::Rosenbrock};
+
+    use super::LBFGSB;
+
+    #[test]
+    fn test_lbfgsb() -> Result<(), Infallible> {
+        let algo = LBFGSB::default();
+        let mut m = Minimizer::new(algo, 2);
+        let problem = Rosenbrock { n: 2 };
+        m.minimize(&problem, &[-2.0, 2.0], &mut ())?;
+        assert!(m.status.converged);
+        assert!(approx_eq!(f64, m.status.fx, 0.0, epsilon = 1e-10));
+        m.minimize(&problem, &[2.0, 2.0], &mut ())?;
+        assert!(m.status.converged);
+        assert!(approx_eq!(f64, m.status.fx, 0.0, epsilon = 1e-10));
+        m.minimize(&problem, &[2.0, -2.0], &mut ())?;
+        assert!(m.status.converged);
+        assert!(approx_eq!(f64, m.status.fx, 0.0, epsilon = 1e-10));
+        m.minimize(&problem, &[-2.0, -2.0], &mut ())?;
+        assert!(m.status.converged);
+        assert!(approx_eq!(f64, m.status.fx, 0.0, epsilon = 1e-10));
+        m.minimize(&problem, &[0.0, 0.0], &mut ())?;
+        assert!(m.status.converged);
+        assert!(approx_eq!(f64, m.status.fx, 0.0, epsilon = 1e-10));
+        m.minimize(&problem, &[1.0, 1.0], &mut ())?;
+        assert!(m.status.converged);
+        assert!(approx_eq!(f64, m.status.fx, 0.0, epsilon = 1e-10));
+        Ok(())
+    }
+}
diff --git a/src/algorithms/nelder_mead.rs b/src/algorithms/nelder_mead.rs
index 0243d00..4c0b0ae 100644
--- a/src/algorithms/nelder_mead.rs
+++ b/src/algorithms/nelder_mead.rs
@@ -712,6 +712,7 @@ where
         bounds: Option<&Vec<Bound<T>>>,
         user_data: &mut U,
     ) -> Result<(), E> {
+        self.status = Status::default();
         self.simplex = self
             .construction_method
             .generate(func, x0, bounds, user_data)?;
@@ -887,8 +888,71 @@ where
             if covariance.is_none() {
                 covariance = hessian.pseudo_inverse(Float::cbrt(T::epsilon())).ok();
             }
-            self.status.cov = covariance
+            self.status.set_cov(covariance);
         }
         Ok(())
     }
 }
+
+#[cfg(test)]
+mod tests {
+    use std::convert::Infallible;
+
+    use float_cmp::approx_eq;
+
+    use crate::{prelude::*, test_functions::Rosenbrock};
+
+    use super::NelderMead;
+
+    #[test]
+    fn test_nelder_mead() -> Result<(), Infallible> {
+        let algo = NelderMead::default();
+        let mut m = Minimizer::new(algo, 2);
+        let problem = Rosenbrock { n: 2 };
+        m.minimize(&problem, &[-2.0, 2.0], &mut ())?;
+        assert!(m.status.converged);
+        assert!(approx_eq!(f64, m.status.fx, 0.0, epsilon = 1e-10));
+        m.minimize(&problem, &[2.0, 2.0], &mut ())?;
+        assert!(m.status.converged);
+        assert!(approx_eq!(f64, m.status.fx, 0.0, epsilon = 1e-10));
+        m.minimize(&problem, &[2.0, -2.0], &mut ())?;
+        assert!(m.status.converged);
+        assert!(approx_eq!(f64, m.status.fx, 0.0, epsilon = 1e-10));
+        m.minimize(&problem, &[-2.0, -2.0], &mut ())?;
+        assert!(m.status.converged);
+        assert!(approx_eq!(f64, m.status.fx, 0.0, epsilon = 1e-10));
+        m.minimize(&problem, &[0.0, 0.0], &mut ())?;
+        assert!(m.status.converged);
+        assert!(approx_eq!(f64, m.status.fx, 0.0, epsilon = 1e-10));
+        m.minimize(&problem, &[1.0, 1.0], &mut ())?;
+        assert!(m.status.converged);
+        assert!(approx_eq!(f64, m.status.fx, 0.0, epsilon = 1e-10));
+        Ok(())
+    }
+
+    #[test]
+    fn test_adaptive_nelder_mead() -> Result<(), Infallible> {
+        let algo = NelderMead::default().with_adaptive(2);
+        let mut m = Minimizer::new(algo, 2);
+        let problem = Rosenbrock { n: 2 };
+        m.minimize(&problem, &[-2.0, 2.0], &mut ())?;
+        assert!(m.status.converged);
+        assert!(approx_eq!(f64, m.status.fx, 0.0, epsilon = 1e-10));
+        m.minimize(&problem, &[2.0, 2.0], &mut ())?;
+        assert!(m.status.converged);
+        assert!(approx_eq!(f64, m.status.fx, 0.0, epsilon = 1e-10));
+        m.minimize(&problem, &[2.0, -2.0], &mut ())?;
+        assert!(m.status.converged);
+        assert!(approx_eq!(f64, m.status.fx, 0.0, epsilon = 1e-10));
+        m.minimize(&problem, &[-2.0, -2.0], &mut ())?;
+        assert!(m.status.converged);
+        assert!(approx_eq!(f64, m.status.fx, 0.0, epsilon = 1e-10));
+        m.minimize(&problem, &[0.0, 0.0], &mut ())?;
+        assert!(m.status.converged);
+        assert!(approx_eq!(f64, m.status.fx, 0.0, epsilon = 1e-10));
+        m.minimize(&problem, &[1.0, 1.0], &mut ())?;
+        assert!(m.status.converged);
+        assert!(approx_eq!(f64, m.status.fx, 0.0, epsilon = 1e-10));
+        Ok(())
+    }
+}
diff --git a/src/lib.rs b/src/lib.rs
index 655e131..8af49ab 100644
--- a/src/lib.rs
+++ b/src/lib.rs
@@ -53,22 +53,31 @@
 //! #             .sum())
 //! #     }
 //! # }
-//!
-//! let problem = Rosenbrock { n: 2 };
-//! let nm = NelderMead::default();
-//! let mut m = Minimizer::new(nm, 2);
-//! let x0 = &[2.0, 2.0];
-//! m.minimize(&problem, x0, &mut ()).unwrap();
-//! println!("{}", m.status);
+//! fn main() -> Result<(), Infallible> {
+//!     let problem = Rosenbrock { n: 2 };
+//!     let nm = NelderMead::default();
+//!     let mut m = Minimizer::new(nm, 2);
+//!     let x0 = &[2.0, 2.0];
+//!     m.minimize(&problem, x0, &mut ())?;
+//!     println!("{}", m.status);
+//!     Ok(())
+//! }
 //! ```
 //!
 //! This should output
 //! ```shell
 //! MSG:       term_f = STDDEV
-//! X:         [0.9999999946231828, 0.9999999884539057]
-//! F(X):      0.00000000000000009170942877687133
-//! N_EVALS:   160
+//! X:         +1.000 ± 0.707
+//!            +1.000 ± 1.416
+//! F(X):      +0.000
+//! N_F_EVALS: 159
+//! N_G_EVALS: 0
 //! CONVERGED: true
+//! COV:       
+//!   ┌             ┐
+//!   │ 0.500 1.000 │
+//!   │ 1.000 2.005 │
+//!   └             ┘
 //! ```
 //!
 //! # Bounds
@@ -500,9 +509,11 @@ impl<T: Scalar> Status<T> {
 }
 impl<T: Scalar + Float> Status<T> {
     /// Sets the covariance matrix and updates parameter errors.
-    pub fn set_cov(&mut self, cov: DMatrix<T>) {
-        self.err = Some(cov.diagonal().map(|v| v.sqrt()));
-        self.cov = Some(cov);
+    pub fn set_cov(&mut self, cov: Option<DMatrix<T>>) {
+        if let Some(cov_mat) = cov {
+            self.err = Some(cov_mat.diagonal().map(|v| v.sqrt()));
+            self.cov = Some(cov_mat);
+        }
     }
 }
 impl<T> Display for Status<T>
@@ -511,19 +522,25 @@ where
 {
     fn fmt(&self, f: &mut std::fmt::Formatter<'_>) -> std::fmt::Result {
         writeln!(f, "MSG:       {}", self.message)?;
-        writeln!(f, "X:         {}", self.x)?;
-        writeln!(f, "F(X):      {}", self.fx)?;
+        write!(f, "X:")?;
+        for i in 0..self.x.len() {
+            if i == 0 {
+                write!(f, "         {:+.3}", self.x[i])?;
+            } else {
+                write!(f, "           {:+.3}", self.x[i])?;
+            }
+            if let Some(e) = &self.err {
+                write!(f, " ± {:.3}", e[i])?;
+            }
+            writeln!(f)?;
+        }
+        writeln!(f, "F(X):      {:+.3}", self.fx)?;
         writeln!(f, "N_F_EVALS: {}", self.n_f_evals)?;
         writeln!(f, "N_G_EVALS: {}", self.n_g_evals)?;
         writeln!(f, "CONVERGED: {}", self.converged)?;
         write!(f, "COV:       ")?;
         match &self.cov {
-            Some(mat) => writeln!(f, "{}", mat),
-            None => writeln!(f, "NOT COMPUTED"),
-        }?;
-        write!(f, "ERR:       ")?;
-        match &self.cov {
-            Some(mat) => writeln!(f, "{}", mat.diagonal().map(|v| v.sqrt())),
+            Some(mat) => writeln!(f, "{:.3}", mat),
             None => writeln!(f, "NOT COMPUTED"),
         }?;
         Ok(())