diff --git a/heat-equation/Cargo.toml b/heat-equation/Cargo.toml
index 8640986..51778c5 100644
--- a/heat-equation/Cargo.toml
+++ b/heat-equation/Cargo.toml
@@ -9,7 +9,7 @@ edition = "2018"
 sbp = { path = "../sbp", features = ["sparse"] }
 ndarray = "0.14.0"
 plotters = { version = "0.3.0", default-features = false, features = ["bitmap_gif", "bitmap_backend", "line_series"] }
-sprs = { version = "0.9.0", default-features = false }
+sprs = { version = "0.10.0", default-features = false }
 
 [dev-dependencies]
 arpack = { git = "https://github.com/mulimoen/arpack-rs", branch = "main" }
diff --git a/heat-equation/src/main.rs b/heat-equation/src/main.rs
index 6b4baba..2ae2cd4 100644
--- a/heat-equation/src/main.rs
+++ b/heat-equation/src/main.rs
@@ -204,10 +204,11 @@ fn dual_dirichlet_sparse(v: ArrayView1<Float>, v0: Float, vn: Float) {
     let rhs = move |fut: &mut Array1<Float>, prev: &Array1<Float>, _t: Float| {
         fut.fill(0.0);
         let prev = prev.as_slice().unwrap();
+        {
+            let fut = fut.as_slice_mut().unwrap();
+            sprs::prod::mul_acc_mat_vec_csr(system.view(), prev, fut);
+        }
         let fut = fut.as_slice_mut().unwrap();
-
-        sprs::prod::mul_acc_mat_vec_csr(system.view(), prev, fut);
-
         sprs::prod::mul_acc_mat_vec_csr(bc.view(), &[v0, vn][..], fut);
     };
 
diff --git a/maxwell/Cargo.toml b/maxwell/Cargo.toml
index 3b108d8..a9fd919 100644
--- a/maxwell/Cargo.toml
+++ b/maxwell/Cargo.toml
@@ -10,7 +10,7 @@ sparse = ["sbp/sparse", "sprs"]
 [dependencies]
 ndarray = "0.14.0"
 sbp = { path = "../sbp" }
-sprs = { version = "0.9.0", optional = true, default-features = false }
+sprs = { version = "0.10.0", optional = true, default-features = false }
 
 [dev-dependencies]
 criterion = "0.3.2"
diff --git a/maxwell/src/sparse.rs b/maxwell/src/sparse.rs
index 42099b6..afc21b3 100644
--- a/maxwell/src/sparse.rs
+++ b/maxwell/src/sparse.rs
@@ -236,11 +236,11 @@ pub fn rhs_matrix(op: &dyn SbpOperator2d, grid: &super::Grid) -> Implicit {
 
 /// RHS with some additional dissipation from the upwind operator
 pub fn rhs_matrix_with_upwind_dissipation(
-    op: &dyn UpwindOperator2d,
+    op: impl UpwindOperator2d + SbpOperator2d,
     grid: &super::Grid,
 ) -> sprs::CsMat<Float> {
-    let rhs = rhs_matrix(op.as_sbp(), grid).rhs;
-    let metrics = grid.metrics(op.as_sbp()).unwrap();
+    let rhs = rhs_matrix(&op, grid).rhs;
+    let metrics = grid.metrics(&op).unwrap();
     let nx = grid.nx();
     let ny = grid.ny();
 
@@ -262,14 +262,14 @@ pub fn rhs_matrix_with_upwind_dissipation(
     };
 
     let diss_x = {
-        let diss_x = UpwindOperator2d::op_xi(op).diss_matrix(nx);
+        let diss_x = UpwindOperator2d::op_xi(&op).diss_matrix(nx);
         let diss_x = kronecker_product(eye(ny).view(), diss_x.view());
         let met = diss(metrics.detj_dxi_dx(), metrics.detj_dxi_dy());
         &met * &kronecker_product(eye(3).view(), diss_x.view())
     };
 
     let diss_y = {
-        let diss_y = UpwindOperator2d::op_eta(op).diss_matrix(ny);
+        let diss_y = UpwindOperator2d::op_eta(&op).diss_matrix(ny);
         let diss_y = kronecker_product(diss_y.view(), eye(nx).view());
         let met = diss(metrics.detj_deta_dx(), metrics.detj_deta_dy());
         &met * &kronecker_product(eye(3).view(), diss_y.view())
diff --git a/sbp/Cargo.toml b/sbp/Cargo.toml
index 72947dd..ac71d9a 100644
--- a/sbp/Cargo.toml
+++ b/sbp/Cargo.toml
@@ -9,7 +9,7 @@ ndarray = { version = "0.14.0", features = ["approx"] }
 approx = "0.4.0"
 packed_simd = { version = "0.3.3", package = "packed_simd_2" }
 rayon = { version = "1.3.0", optional = true }
-sprs = { version = "0.9.0", optional = true, default-features = false }
+sprs = { version = "0.10.0", optional = true, default-features = false }
 serde = { version = "1.0.115", optional = true, default-features = false, features = ["derive"] }
 num-traits = "0.2.14"
 
diff --git a/sbp/src/operators/traditional4.rs b/sbp/src/operators/traditional4.rs
index 2bb2ed6..9022635 100644
--- a/sbp/src/operators/traditional4.rs
+++ b/sbp/src/operators/traditional4.rs
@@ -111,12 +111,13 @@ impl super::SbpOperator1d2 for SBP4 {
             // d_x => -d_x when reversed
             d1.iter_mut().for_each(|v| *v *= -1.0);
             // Indices are now in reverse order
-            sprs::CsMat::new(
+            sprs::CsMat::new_from_unsorted(
                 (1, n),
                 vec![0, d1.len()],
                 (0..n).rev().take(d1.len()).collect(),
                 d1,
             )
+            .unwrap()
         }
     }
 }