diff --git a/qmat/lagrange.py b/qmat/lagrange.py
index 63300a0..b34a198 100644
--- a/qmat/lagrange.py
+++ b/qmat/lagrange.py
@@ -530,3 +530,51 @@ def getDerivationMatrix(self, order=1, duplicates=True):
             return D2
         else:
             return D1, D2
+
+
+def getSparseInterpolationMatrix(inPoints, outPoints, order):
+    """
+    Get a sparse interpolation matrix from `inPoints` to `outPoints` of order
+    `order` using barycentric Lagrange interpolation.
+
+    The matrix will have `order` entries per line, and tends to be banded when
+    both `inPoints` and `outPoints` are equispaced and cover the same interval.
+
+    Parameters
+    ----------
+        inPoints : np.1darray
+            The points you want to interpolate from
+        outPoints : np.1darray
+            The points you want to interpolate to
+        order : int
+            Order of the interpolation
+
+    Returns
+    -------
+    A : scipy.sparse.csc_matrix(len(outPoints), len(inPoints))
+        Sparse interpolation matrix
+    """
+    import scipy.sparse as sp
+
+    assert order <= len(inPoints), f'Cannot interpolate {len(inPoints)} to order {order}! Please reduce order'
+
+    A = sp.lil_matrix((len(outPoints), len(inPoints)))
+    lastInterpolationLine = None
+    lastClosestPoints = None
+
+    for i in range(len(outPoints)):
+        closestPointsIdx = np.sort(np.argsort(np.abs(inPoints - outPoints[i]))[:order])
+        closestPoints = inPoints[closestPointsIdx] - outPoints[i]
+
+        if lastClosestPoints is not None and np.allclose(closestPoints, lastClosestPoints):
+            interpolationLine = lastInterpolationLine
+        else:
+            interpolator = LagrangeApproximation(points = closestPoints)
+            interpolationLine = interpolator.getInterpolationMatrix([0])[0]
+
+            lastInterpolationLine = interpolationLine
+            lastClosestPoints = closestPoints
+
+        A[i, closestPointsIdx] = interpolationLine
+    
+    return A.tocsc()
diff --git a/tests/test_2_lagrange.py b/tests/test_2_lagrange.py
index b0e0674..dc3e53a 100644
--- a/tests/test_2_lagrange.py
+++ b/tests/test_2_lagrange.py
@@ -166,4 +166,34 @@ def testDuplicates(nPoints, nCopy, duplicates):
 
     assert np.allclose(D2[:, approx._zerIdx], 0), "[D2] zero indices have non-zero values"
     assert np.allclose(D2[:, approx._nnzIdx], D2_ref), "[D2] nonzero values different from reference"
-    assert np.allclose(D2_noDuplicates, D2_ref), "[D2] no duplicates values different from reference"
\ No newline at end of file
+    assert np.allclose(D2_noDuplicates, D2_ref), "[D2] no duplicates values different from reference"
+
+
+@pytest.mark.parametrize('inPoints', [np.linspace(0, 1, 15), np.linspace(0.3, 2.7, 13), np.linspace(0, 10, 4)])
+@pytest.mark.parametrize('outPoints', [np.linspace(0, 1, 9), np.linspace(0.5, 2.2, 4), np.linspace(0.3, 2.7, 13)])
+@pytest.mark.parametrize('order', [1, 2, 3, 4])
+def testSparseInterpolation(inPoints, outPoints, order):
+    from qmat.lagrange import getSparseInterpolationMatrix
+    import scipy.sparse as sp
+
+    np.random.seed(47)
+
+    interpolationMatrix = getSparseInterpolationMatrix(inPoints, outPoints, order)
+    assert isinstance(interpolationMatrix, sp.csc_matrix)
+
+    polyCoeffs = np.random.randn(order)
+    inPolynomial = np.polyval(polyCoeffs,inPoints)
+    interpolated = interpolationMatrix @ inPolynomial
+    error = np.linalg.norm(np.polyval(polyCoeffs, outPoints)- interpolated)
+    assert error < 1e-11, f'Interpolation of order {order} polynomial is not exact with error {error:.2e}'
+
+    testInexactInterpolation = True
+    if len(inPoints) == len(outPoints):
+        if np.allclose(inPoints, outPoints):
+            testInexactInterpolation = False
+
+    if testInexactInterpolation:
+        polyCoeffs = np.random.randn(order+1)
+        inPolynomial = np.polyval(polyCoeffs,inPoints)
+        interpolated = interpolationMatrix @ inPolynomial
+        assert not np.allclose(np.polyval(polyCoeffs, outPoints), interpolated), f'Interpolation of order {order+1} polynomial is unexpectedly exact'