diff --git a/pyproject.toml b/pyproject.toml
index afe3feb4..98af72ab 100644
--- a/pyproject.toml
+++ b/pyproject.toml
@@ -35,7 +35,6 @@ classifiers = [
 ]
 requires-python = ">=3.10"
 dependencies = [
-    "findiff>=0.9.2",
     "h5py>=3.8.0",
     "numpy>=1.21.5",
     "scipy>=1.9.0",
diff --git a/src/WallGo/boltzmann.py b/src/WallGo/boltzmann.py
index 5e311a5d..267e1eca 100644
--- a/src/WallGo/boltzmann.py
+++ b/src/WallGo/boltzmann.py
@@ -8,7 +8,6 @@
 import logging
 from enum import Enum, auto
 import numpy as np
-import findiff  # finite difference methods
 from .containers import BoltzmannBackground, BoltzmannDeltas
 from .grid import Grid
 from .polynomial import Polynomial, SpectralConvergenceInfo
@@ -16,6 +15,7 @@
 from .collisionArray import CollisionArray
 from .results import BoltzmannResults
 from .exceptions import CollisionLoadError
+from .helpers import finiteDiffMatrix
 
 if typing.TYPE_CHECKING:
     import importlib
@@ -704,10 +704,8 @@ def buildLinearEquations(
             intertwinerRpMat = np.identity(self.grid.N - 1)
             # derivative matrices
             chiFull, rzFull, _ = self.grid.getCompactCoordinates(endpoints=True)
-            derivOperatorChi = findiff.FinDiff((0, chiFull, 1), acc=2)
-            derivMatrixChi = derivOperatorChi.matrix((self.grid.M + 1,))
-            derivOperatorRz = findiff.FinDiff((0, rzFull, 1), acc=2)
-            derivMatrixRz = derivOperatorRz.matrix((self.grid.N + 1,))
+            derivMatrixChi = finiteDiffMatrix(chiFull)
+            derivMatrixRz = finiteDiffMatrix(rzFull)
             # spatial derivatives of profiles, endpoints used for taking
             # derivatives but then dropped as deltaF fixed at 0 at endpoints
             dTemperaturedChi = (derivMatrixChi @ temperatureFull)[
diff --git a/src/WallGo/helpers.py b/src/WallGo/helpers.py
index eeef79a5..b6e1eee8 100644
--- a/src/WallGo/helpers.py
+++ b/src/WallGo/helpers.py
@@ -425,6 +425,34 @@ def hessian(
     fEvaluation = f(pos, *args).reshape(shape)
     return np.asarray(np.sum(coeff * fEvaluation, axis=-3))
 
+def finiteDiffMatrix(grid: np.ndarray) -> np.ndarray:
+    """
+    Computes the finite difference matrix corresponding to the array grid.
+    grid must contain the endpoints.
+    """
+    diffMatrix = np.zeros((grid.size, grid.size))
+    dx1 = -(grid[1:-1] - grid[:-2])
+    dx2 = grid[2:] - grid[1:-1]
+
+    # Below the diagonal
+    diffMatrix[1:-1,:-2] += np.diag(dx2/(dx1*(dx2-dx1)))
+    # Diagonal
+    diffMatrix[1:-1,1:-1] += np.diag(-(dx1+dx2)/(dx1*dx2))
+    # Above the diagonal
+    diffMatrix[1:-1,2:] += np.diag(dx1/(dx2*(dx1-dx2)))
+
+    # Left boundary
+    diffMatrix[0,0] = ((2*grid[0] - grid[1] - grid[2])/((grid[0] - grid[1])*(grid[0] - grid[2])))
+    diffMatrix[0,1] = ((-grid[0] + grid[2])/((grid[0] - grid[1])*(grid[1] - grid[2])))
+    diffMatrix[0,2] = ((-grid[0] + grid[1])/((grid[0] - grid[2])*(-grid[1] + grid[2])))
+
+    # Right boundary
+    diffMatrix[-1,-1] = ((2*grid[-1] - grid[-2] - grid[-3])/((grid[-1] - grid[-2])*(grid[-1] - grid[-3])))
+    diffMatrix[-1,-2] = ((-grid[-1] + grid[-3])/((grid[-1] - grid[-2])*(grid[-2] - grid[-3])))
+    diffMatrix[-1,-3] = ((-grid[-1] + grid[-2])/((grid[-1] - grid[-3])*(-grid[-2] + grid[-3])))
+
+    return diffMatrix
+
 
 def gammaSq(v: float) -> float:
     r"""