diff --git a/test/usecases/solve/finite_difference.mod b/test/usecases/solve/finite_difference.mod
new file mode 100644
index 000000000..2c0f94e86
--- /dev/null
+++ b/test/usecases/solve/finite_difference.mod
@@ -0,0 +1,30 @@
+NEURON {
+    SUFFIX finite_difference
+    GLOBAL a
+    THREADSAFE
+}
+
+ASSIGNED {
+  a
+}
+
+STATE {
+  x
+}
+
+INITIAL {
+    x = 42.0
+    a = 0.1
+}
+
+BREAKPOINT {
+    SOLVE dX METHOD derivimplicit
+}
+
+DERIVATIVE dX {
+    x' = -f(x)
+}
+
+FUNCTION f(x) {
+    f = a*x
+}
diff --git a/test/usecases/solve/test_finite_difference.py b/test/usecases/solve/test_finite_difference.py
new file mode 100644
index 000000000..ba204ff2f
--- /dev/null
+++ b/test/usecases/solve/test_finite_difference.py
@@ -0,0 +1,30 @@
+import numpy as np
+from neuron import h, gui
+from neuron.units import ms
+
+
+def test_finite_difference():
+    nseg = 1
+
+    s = h.Section()
+    s.insert("finite_difference")
+    s.nseg = nseg
+
+    x_hoc = h.Vector().record(getattr(s(0.5), f"_ref_x_finite_difference"))
+    t_hoc = h.Vector().record(h._ref_t)
+
+    h.stdinit()
+    h.dt = 0.001
+    h.tstop = 5.0 * ms
+    h.run()
+
+    x = np.array(x_hoc.as_numpy())
+    t = np.array(t_hoc.as_numpy())
+
+    a = h.a_finite_difference
+    x_exact = 42.0 * np.exp(-a * t)
+    np.testing.assert_allclose(x, x_exact, rtol=1e-4)
+
+
+if __name__ == "__main__":
+    test_finite_difference()