diff --git a/examples/solveh.rs b/examples/solveh.rs
index 506a1aa7..26214d6f 100644
--- a/examples/solveh.rs
+++ b/examples/solveh.rs
@@ -10,8 +10,7 @@ fn solve() -> Result<(), error::LinalgError> {
     let a: Array2<c64> = random_hermite(3); // complex Hermite positive definite matrix
     let b: Array1<c64> = random(3);
     println!("b = {:?}", &b);
-    let f = a.factorizeh()?; // DK factorize
-    let x = f.solveh(b)?;
+    let x = a.solveh(&b)?;
     println!("Ax = {:?}", a.dot(&x));;
     Ok(())
 }
@@ -23,7 +22,7 @@ fn factorize() -> Result<(), error::LinalgError> {
     // once factorized, you can use it several times:
     for _ in 0..10 {
         let b: Array1<f64> = random(3);
-        let _x = f.solveh(b)?;
+        let _x = f.solveh_into(b)?;
     }
     Ok(())
 }
diff --git a/src/solveh.rs b/src/solveh.rs
index e72a3f19..62c3728f 100644
--- a/src/solveh.rs
+++ b/src/solveh.rs
@@ -49,6 +49,21 @@ where
     }
 }
 
+impl<A, S> SolveH<A> for ArrayBase<S, Ix2>
+where
+    A: Scalar,
+    S: Data<Elem = A>,
+{
+    fn solveh_mut<'a, Sb>(&self, mut rhs: &'a mut ArrayBase<Sb, Ix1>) -> Result<&'a mut ArrayBase<Sb, Ix1>>
+    where
+        Sb: DataMut<Elem = A>,
+    {
+        let f = self.factorizeh()?;
+        f.solveh_mut(rhs)
+    }
+}
+
+
 impl<A, S> FactorizedH<S>
 where
     A: Scalar,