diff --git a/sbp/src/operators.rs b/sbp/src/operators.rs
index 1457429..83090ab 100644
--- a/sbp/src/operators.rs
+++ b/sbp/src/operators.rs
@@ -78,3 +78,47 @@ mod traditional4;
 pub use traditional4::SBP4;
 mod traditional8;
 pub use traditional8::SBP8;
+
+#[cfg(test)]
+pub(crate) mod testing {
+    use super::*;
+    use ndarray::prelude::*;
+    pub(crate) fn grid_eval<F: Fn(Float, Float) -> Float>(
+        n: (usize, usize),
+        f: F,
+    ) -> Array2<Float> {
+        let nx = n.1;
+        let dx = 1.0 / (nx - 1) as Float;
+        let ny = n.0;
+        let dy = 1.0 / (ny - 1) as Float;
+        Array2::from_shape_fn(n, |(j, i)| {
+            let x = dx * i as Float;
+            let y = dy * j as Float;
+            f(x, y)
+        })
+    }
+
+    pub(crate) fn check_operator_on<SBP, F, FX, FY>(
+        n: (usize, usize),
+        f: F,
+        dfdx: FX,
+        dfdy: FY,
+        eps: Float,
+    ) where
+        SBP: SbpOperator,
+        F: Fn(Float, Float) -> Float,
+        FX: Fn(Float, Float) -> Float,
+        FY: Fn(Float, Float) -> Float,
+    {
+        let mut y = Array2::zeros(n);
+        let x = grid_eval(n, f);
+
+        y.fill(0.0);
+        SBP::diffxi(x.view(), y.view_mut());
+        approx::assert_abs_diff_eq!(&y, &grid_eval(n, dfdx), epsilon = eps);
+
+        y.fill(0.0);
+        SBP::diffeta(x.view(), y.view_mut());
+        approx::assert_abs_diff_eq!(&y, &grid_eval(n, dfdy), epsilon = eps);
+    }
+}
diff --git a/sbp/src/operators/traditional4.rs b/sbp/src/operators/traditional4.rs
index f1b546b..91706d4 100644
--- a/sbp/src/operators/traditional4.rs
+++ b/sbp/src/operators/traditional4.rs
@@ -45,3 +45,27 @@ impl SbpOperator for SBP4 {
         Self::HBLOCK
     }
 }
+
+#[test]
+fn test_trad4() {
+    use super::testing::*;
+    use super::*;
+    let nx = 20;
+    let ny = 13;
+
+    check_operator_on::<SBP4, _, _, _>((ny, nx), |x, y| x + 2.0 * y, |_, _| 1.0, |_, _| 2.0, 1e-4);
+    check_operator_on::<SBP4, _, _, _>(
+        (ny, nx),
+        |x, y| x * x + 2.0 * x * y + 3.0 * y * y,
+        |x, y| 2.0 * x + 2.0 * y,
+        |x, y| 2.0 * x + 6.0 * y,
+        1e-3,
+    );
+    check_operator_on::<SBP4, _, _, _>(
+        (ny, nx),
+        |x, y| x.powi(3) + 2.0 * x.powi(2) * y + 3.0 * x * y.powi(2) + 4.0 * y.powi(3),
+        |x, y| 3.0 * x.powi(2) + 4.0 * x * y + 3.0 * y.powi(2),
+        |x, y| 2.0 * x.powi(2) + 6.0 * x * y + 12.0 * y.powi(2),
+        1e-1,
+    );
+}
diff --git a/sbp/src/operators/traditional8.rs b/sbp/src/operators/traditional8.rs
index fb390cc..856e367 100644
--- a/sbp/src/operators/traditional8.rs
+++ b/sbp/src/operators/traditional8.rs
@@ -49,3 +49,61 @@ impl SbpOperator for SBP8 {
         Self::HBLOCK
     }
 }
+
+#[test]
+fn test_trad8() {
+    use super::testing::*;
+    let nx = 32;
+    let ny = 16;
+
+    // Order one polynomial
+    check_operator_on::<SBP8, _, _, _>(
+        (ny, nx),
+        |x, y| x + 2.0 * y,
+        |_x, _y| 1.0,
+        |_x, _y| 2.0,
+        1e-4,
+    );
+
+    // Order two polynomial
+    check_operator_on::<SBP8, _, _, _>(
+        (ny, nx),
+        |x, y| x * x + 0.5 * y * y,
+        |x, _y| 2.0 * x,
+        |_x, y| y,
+        1e-4,
+    );
+    check_operator_on::<SBP8, _, _, _>((ny, nx), |x, y| x * y, |_x, y| y, |x, _y| x, 1e-4);
+
+    // Order three polynomials
+    check_operator_on::<SBP8, _, _, _>(
+        (ny, nx),
+        |x, y| x * x * x + y * y * y / 6.0,
+        |x, _y| 3.0 * x * x,
+        |_x, y| y * y / 2.0,
+        1e-4,
+    );
+    check_operator_on::<SBP8, _, _, _>(
+        (ny, nx),
+        |x, y| x * x * y + x * y * y / 2.0,
+        |x, y| 2.0 * x * y + y * y / 2.0,
+        |x, y| x * x + x * y,
+        1e-4,
+    );
+    check_operator_on::<SBP8, _, _, _>(
+        (ny, nx),
+        |x, y| x.powi(3) + 2.0 * x.powi(2) * y + 3.0 * x * y.powi(2) + 4.0 * y.powi(3),
+        |x, y| 3.0 * x.powi(2) + 4.0 * x * y + 3.0 * y.powi(2),
+        |x, y| 2.0 * x.powi(2) + 6.0 * x * y + 12.0 * y.powi(2),
+        1e-4,
+    );
+
+    // Order four polynomials
+    check_operator_on::<SBP8, _, _, _>(
+        (ny, nx),
+        |x, y| x.powi(4) + x.powi(3) * y + x.powi(2) * y.powi(2) + x * y.powi(3) + y.powi(4),
+        |x, y| 4.0 * x.powi(3) + 3.0 * x.powi(2) * y + 2.0 * x * y.powi(2) + y.powi(3),
+        |x, y| x.powi(3) + 2.0 * x.powi(2) * y + 3.0 * x * y.powi(2) + 4.0 * y.powi(3),
+        1e-4,
+    );
+}
diff --git a/sbp/src/operators/upwind4.rs b/sbp/src/operators/upwind4.rs
index 44d3b12..0de501b 100644
--- a/sbp/src/operators/upwind4.rs
+++ b/sbp/src/operators/upwind4.rs
@@ -451,3 +451,40 @@ impl UpwindOperator for Upwind4 {
         Self::dissxi(prev.reversed_axes(), fut.reversed_axes());
     }
 }
+
+#[test]
+fn upwind4_test2() {
+    use super::testing::*;
+    use super::*;
+    let nx = 32;
+    let ny = 16;
+
+    check_operator_on::<Upwind4, _, _, _>(
+        (ny, nx),
+        |x, y| x + 2.0 * y,
+        |_, _| 1.0,
+        |_, _| 2.0,
+        1e-4,
+    );
+    check_operator_on::<Upwind4, _, _, _>(
+        (ny, nx),
+        |x, y| x * x + 2.0 * x * y + 3.0 * y * y,
+        |x, y| 2.0 * x + 2.0 * y,
+        |x, y| 2.0 * x + 6.0 * y,
+        1e-3,
+    );
+    check_operator_on::<Upwind4, _, _, _>(
+        (ny, nx),
+        |x, y| x.powi(3) + 2.0 * x.powi(2) * y + 3.0 * x * y.powi(2) + 4.0 * y.powi(3),
+        |x, y| 3.0 * x.powi(2) + 4.0 * x * y + 3.0 * y.powi(2),
+        |x, y| 2.0 * x.powi(2) + 6.0 * x * y + 12.0 * y.powi(2),
+        1e-1,
+    );
+    check_operator_on::<Upwind4, _, _, _>(
+        (32, 32),
+        |x, y| x.powi(3) + 2.0 * x.powi(2) * y + 3.0 * x * y.powi(2) + 4.0 * y.powi(3),
+        |x, y| 3.0 * x.powi(2) + 4.0 * x * y + 3.0 * y.powi(2),
+        |x, y| 2.0 * x.powi(2) + 6.0 * x * y + 12.0 * y.powi(2),
+        1e-1,
+    );
+}
diff --git a/sbp/src/operators/upwind9.rs b/sbp/src/operators/upwind9.rs
index d5cb69b..3b11b99 100644
--- a/sbp/src/operators/upwind9.rs
+++ b/sbp/src/operators/upwind9.rs
@@ -89,3 +89,61 @@ impl UpwindOperator for Upwind9 {
         Self::dissxi(prev.reversed_axes(), fut.reversed_axes());
     }
 }
+
+#[test]
+fn test_upwind9() {
+    use super::testing::*;
+    let nx = 32;
+    let ny = 16;
+
+    // Order one polynomial
+    check_operator_on::<Upwind9, _, _, _>(
+        (ny, nx),
+        |x, y| x + 2.0 * y,
+        |_x, _y| 1.0,
+        |_x, _y| 2.0,
+        1e-4,
+    );
+
+    // Order two polynomial
+    check_operator_on::<Upwind9, _, _, _>(
+        (ny, nx),
+        |x, y| x * x + 0.5 * y * y,
+        |x, _y| 2.0 * x,
+        |_x, y| y,
+        1e-4,
+    );
+    check_operator_on::<Upwind9, _, _, _>((ny, nx), |x, y| x * y, |_x, y| y, |x, _y| x, 1e-4);
+
+    // Order three polynomials
+    check_operator_on::<Upwind9, _, _, _>(
+        (ny, nx),
+        |x, y| x * x * x + y * y * y / 6.0,
+        |x, _y| 3.0 * x * x,
+        |_x, y| y * y / 2.0,
+        1e-4,
+    );
+    check_operator_on::<Upwind9, _, _, _>(
+        (ny, nx),
+        |x, y| x * x * y + x * y * y / 2.0,
+        |x, y| 2.0 * x * y + y * y / 2.0,
+        |x, y| x * x + x * y,
+        1e-4,
+    );
+    check_operator_on::<Upwind9, _, _, _>(
+        (ny, nx),
+        |x, y| x.powi(3) + 2.0 * x.powi(2) * y + 3.0 * x * y.powi(2) + 4.0 * y.powi(3),
+        |x, y| 3.0 * x.powi(2) + 4.0 * x * y + 3.0 * y.powi(2),
+        |x, y| 2.0 * x.powi(2) + 6.0 * x * y + 12.0 * y.powi(2),
+        1e-4,
+    );
+
+    // Order four polynomials
+    check_operator_on::<Upwind9, _, _, _>(
+        (ny, nx),
+        |x, y| x.powi(4) + x.powi(3) * y + x.powi(2) * y.powi(2) + x * y.powi(3) + y.powi(4),
+        |x, y| 4.0 * x.powi(3) + 3.0 * x.powi(2) * y + 2.0 * x * y.powi(2) + y.powi(3),
+        |x, y| x.powi(3) + 2.0 * x.powi(2) * y + 3.0 * x * y.powi(2) + 4.0 * y.powi(3),
+        1e-4,
+    );
+}