diff --git a/src/pycaputo/differentiation/caputo.py b/src/pycaputo/differentiation/caputo.py
index 5740d71..a9144e7 100644
--- a/src/pycaputo/differentiation/caputo.py
+++ b/src/pycaputo/differentiation/caputo.py
@@ -64,25 +64,37 @@ def __post_init__(self) -> None:
                 )
 
 
-@quadrature_weights.register(L1)
-def _quadrature_weights_caputo_l1(m: L1, p: Points, n: int) -> Array:
-    if n < 0:
-        n = p.size + n
+def _caputo_piecewise_constant_integral(p: Points, n: int, alpha: float) -> Array:
+    r"""Computes the integrals
 
-    if not 0 <= n <= p.size:
-        raise IndexError(f"Index 'n' out of range: 0 <= {n} < {p.size}")
+    .. math::
+
+        a_{nk} = \int_{x_k}^{x_{k + 1}} \frac{1}{(x_n - s)^\alpha} ds
+
+    for :math:`0 \le k \le n - 1`.
+    """
 
-    alpha = m.alpha
     x, dx = p.x[:n], p.dx[: n - 1]
     xn = x[-1]
 
-    w = np.empty_like(x)
-    w[n - 1] = 0.0
-    w[: n - 1] = (
+    return (
         ((xn - x[:-1]) ** (1 - alpha) - (xn - x[1:]) ** (1 - alpha))
         / dx
         / math.gamma(2 - alpha)
     )
+
+
+@quadrature_weights.register(L1)
+def _quadrature_weights_caputo_l1(m: L1, p: Points, n: int) -> Array:
+    if n < 0:
+        n = p.size + n
+
+    if not 0 <= n <= p.size:
+        raise IndexError(f"Index 'n' out of range: 0 <= {n} < {p.size}")
+
+    w = np.empty(n, dtype=p.x.dtype)
+    w[n - 1] = 0.0
+    w[: n - 1] = _caputo_piecewise_constant_integral(p, n, m.alpha)
     w[1:n] = w[: n - 1] - w[1:n]
     w[0] = -w[0]
 
@@ -97,27 +109,22 @@ def _diffs_caputo_l1(m: L1, f: ArrayOrScalarFunction, p: Points, n: int) -> Scal
     if not 0 <= n <= p.size:
         raise IndexError(f"Index 'n' out of range: 0 <= {n} < {p.size}")
 
-    if n == 1:
+    if n == 0:
         return np.array([np.nan])
 
     w = quadrature_weights(m, p, n + 1)
-    fx = f(p.x[: n + 1]) if callable(f) else f
+    fx = f(p.x[: n + 1]) if callable(f) else f[:n + 1]
 
     return np.sum(w * fx)  # type: ignore[no-any-return]
 
 
 @diff.register(L1)
 def _diff_caputo_l1(m: L1, f: ArrayOrScalarFunction, p: Points) -> Array:
-    if callable(f):
-        fx = f(p.x)
-    elif isinstance(f, np.ndarray):
-        if f.shape[0] != p.size:
-            raise ValueError(
-                f"Shape of 'f' does match points: got {f.shape} expected{p.shape}"
-            )
-        fx = f
-    else:
-        raise TypeError(f"Unknown type for 'f': {type(f)}")
+    fx = f(p.x) if callable(f) else f
+    if fx.shape[0] != p.size:
+        raise ValueError(
+            f"Shape of 'f' does match points: got {f.shape} expected {p.shape}"
+        )
 
     df = np.empty_like(fx)
     df[0] = np.nan
@@ -145,6 +152,10 @@ class ModifiedL1(CaputoMethod):
     This method is defined in Section 4.1.1 (III) from [Li2020]_. Note that this
     method evaluates the fractional derivative at the midpoints
     :math:`(x_i + x_{i - 1}) / 2` for any given input grid.
+
+    As noted in [Li2020]_, this method has the same order of convergence as the
+    standard L1 method. However the weights can be used to construct a
+    Crank-Nicolson type method.
     """
 
     if __debug__:
@@ -225,6 +236,90 @@ def __post_init__(self) -> None:
                 )
 
 
+def _caputo_d2_coefficients(x: Array) -> Array:
+    x0, x1, x2, x3 = x[:4]
+    a0 = 2.0 * (3.0 * x0 - x1 - x2 - x3) / ((x0 - x1) * (x0 - x2) * (x0 - x3))
+    b0 = 2.0 * (x3 - 2.0 * x0 + x2) / ((x0 - x1) * (x1 - x2) * (x1 - x3))
+    c0 = 2.0 * (x3 - 2.0 * x0 + x1) / ((x0 - x2) * (x2 - x1) * (x2 - x3))
+    d0 = 2.0 * (x2 - 2.0 * x0 + x1) / ((x0 - x3) * (x3 - x1) * (x3 - x2))
+
+    return np.array([a0, b0, c0, d0])
+
+
+@quadrature_weights.register(L2)
+def _quadrature_weights_caputo_l2(m: L2, p: Points, n: int) -> Array:
+    if n < 0:
+        n = p.size + n
+
+    if not 0 <= n <= p.size:
+        raise IndexError(f"Index 'n' out of range: 0 <= {n} < {p.size}")
+
+    dx = p.dx[:n - 1]
+    dxm = (dx[1:] + dx[:-1]) / 2.0
+
+    # get finite difference coefficients for center stencil
+    a = 1.0 / (dx[:-2] * dxm)
+    b = 1.0 / (dx[1:-1] * dxm)
+
+    # get finite difference coefficients for boundary stencils
+    a_l, b_l, c_l, d_l = _caputo_d2_coefficients(p.x)
+    d_r, c_r, b_r, a_r = _caputo_d2_coefficients(p.x[::-1])
+
+    wi = _caputo_piecewise_constant_integral(p, n, m.alpha)
+
+    w = np.empty_like(p.x[:n])
+    # center stencil
+    w[2:n - 1] = a * wi[2:] - (a + b) * w[1:-1] + b * w[:-2]
+    # add left boundary stencil
+    w[0] = a_l + a[1] * wi[1]
+    w[1] = b_l + a[2] * wi[2] - (a[1] + b[1]) * wi[1]
+    w[2] += c_l
+    w[3] += d_l
+    # add right-boundary stencil
+    if n == p.size:
+        w[-4] += d_r
+        w[-3] += c_r
+        w[-2] = b_r + a[-2] * wi[-2] - (a[-1] + b[-1]) * wi[-1]
+        w[-1] = a_r + a[-1] * wi[-1]
+
+    return w
+
+
+@diffs.register(L2)
+def _diffs_caputo_l2(m: L2, f: ArrayOrScalarFunction, p: Points, n: int) -> Array:
+    if n < 0:
+        n = p.size + n
+
+    if not 0 <= n <= p.size:
+        raise IndexError(f"Index 'n' out of range: 0 <= {n} < {p.size}")
+
+    if n == 0:
+        return np.array([np.nan])
+
+    w = quadrature_weights(m, p, n + 1)
+    fx = f(p.x[: n + 1]) if callable(f) else f[:n + 1]
+
+    return np.sum(w * fx)  # type: ignore[no-any-return]
+
+
+@diff.register(L2)
+def _diff_caputo_l2(m: L2, f: ArrayOrScalarFunction, p: Points) -> Array:
+    fx = f(p.x) if callable(f) else f
+    if fx.shape[0] != p.size:
+        raise ValueError(
+            f"Shape of 'f' does match points: got {f.shape} expected {p.shape}"
+        )
+
+    df = np.empty_like(fx)
+    df[0] = np.nan
+
+    # FIXME: in the uniform case, we can also do an FFT, but we need different
+    # weights for that, so we leave it like this for now
+    for n in range(1, df.size):
+        w = quadrature_weights(m, p, n + 1)
+        df[n] = np.sum(w * fx[: n + 1])
+
+
 def _weights_l2(alpha: float, i: int | Array, k: int | Array) -> Array:
     return np.array((i - k) ** (2 - alpha) - (i - k - 1) ** (2 - alpha))
 
@@ -377,7 +472,7 @@ def _diff_jacobi(m: SpectralJacobi, f: ArrayOrScalarFunction, p: Points) -> Arra
 # }}}
 
 
-# {{{
+# {{{ DiffusiveCaputoMethod
 
 
 @dataclass(frozen=True)
diff --git a/tests/test_diff_caputo.py b/tests/test_diff_caputo.py
index 881fb23..3a6881c 100644
--- a/tests/test_diff_caputo.py
+++ b/tests/test_diff_caputo.py
@@ -65,6 +65,7 @@ def df_test(x: Array, *, alpha: float, mu: float = 3.5) -> Array:
         ("L2", "uniform"),
         ("L2C", "uniform"),
         ("ModifiedL1", "midpoints"),
+        ("ModifiedL1", "stretch"),
     ],
 )
 @pytest.mark.parametrize("alpha", [0.1, 0.25, 0.5, 0.75, 0.9])
@@ -114,6 +115,11 @@ def test_caputo_lmethods(
 
     for n in [16, 32, 64, 128, 256, 512, 768, 1024]:
         p = make_points_from_name(grid_type, n, a=0.0, b=0.5)
+        if name == "ModifiedL1" and grid_type != "midpoints":
+            from pycaputo.grid import make_midpoints_from
+
+            p = make_midpoints_from(p)
+
         df_num = diff(meth, f_test, p)
         df_ref = df_test(p.x, alpha=alpha)