diff --git a/optimism/Interpolants.py b/optimism/Interpolants.py
index 76f1dee..63a424c 100644
--- a/optimism/Interpolants.py
+++ b/optimism/Interpolants.py
@@ -81,7 +81,7 @@ def make_parent_element_1d(degree):
 def get_lobatto_nodes_1d(degree):
     p = onp.polynomial.Legendre.basis(degree, domain=[0.0, 1.0])
     dp = p.deriv()
-    xInterior = dp.roots()
+    xInterior = np.real(dp.roots())
     xn = np.hstack((np.array([0.0]), xInterior, np.array([1.0])))
     return xn
 
diff --git a/optimism/SparseCholesky.py b/optimism/SparseCholesky.py
index e6e4358..44206ee 100644
--- a/optimism/SparseCholesky.py
+++ b/optimism/SparseCholesky.py
@@ -3,7 +3,7 @@
 
 import numpy as onp
 
-from sksparse.cholmod import analyze, cholesky
+from sksparse import cholmod
 from sksparse.cholmod import CholmodNotPositiveDefiniteError as NotPosDefError
 
 from optimism.JaxConfig import *
@@ -21,15 +21,15 @@ def factorize(self):
         # we can improve this later if we are inclined
         assert isspmatrix_csc(self.A),  \
             "Preconditioner matrix is not in a valid sparse format"
-        self.Precond = analyze(self.A, mode='supernodal',
-                               ordering_method='nesdis')
+        self.Precond = cholmod.CholeskyFactor(self.A, sym_kind="sym", supernodal_mode="supernodal",
+                                              order="metis")
 
         attempt = 0
         maxAttempts = 10
         while attempt < maxAttempts:
             try:
                 print('Factorizing preconditioner')
-                self.Precond.cholesky_inplace(self.A)
+                self.Precond.factorize(self.A)
             except NotPosDefError:
                 attempt += 1
                 print('Cholesky failed, assembling preconditioner', attempt)
@@ -41,7 +41,7 @@ def factorize(self):
         if attempt == maxAttempts:
             print("Cholesky failed too many times, using identity preconditioner")
             self.A = identity(self.A.shape[0], format='csc')
-            self.Precond.cholesky_inplace(self.A)
+            self.Precond.factorize(self.A)
 
             
     def update(self, new_stiffness_func):
@@ -50,9 +50,7 @@ def update(self, new_stiffness_func):
 
         
     def apply(self, b):
-        if type(b) == type(np.array([])):
-            b = onp.array(b, copy=False)
-        return self.Precond(b)
+        return np.asarray(self.Precond.solve(onp.array(b, copy=True)))
 
         
     def apply_transpose(self, b):
@@ -67,16 +65,6 @@ def multiply_by_transpose(self, x):
         return self.A.T.dot(x)
 
 
-    def check_stability(self, x, p):
-        A = self.stiffness_func(x, p)
-        try:
-            self.Precond.cholesky(A)
-            print("Jacobian is stable.")
-        except NotPosDefError as e:
-            print(e)
-            print("Jacobian is unstable.")
-
-
     def get_diagonal_stiffness(self):
         return self.A.diagonal()
 
diff --git a/optimism/test/test_TensorMath.py b/optimism/test/test_TensorMath.py
index eb69527..c365767 100644
--- a/optimism/test/test_TensorMath.py
+++ b/optimism/test/test_TensorMath.py
@@ -265,7 +265,7 @@ def test_right_polar_decomp(self):
         # U is symmetric
         self.assertArrayNear(U, TensorMath.sym(U), 14)
         # RU = F
-        self.assertArrayNear(R@U, F, 14)
+        self.assertArrayNear(R@U, F, 13)
 
 
 if __name__ == '__main__':