diff --git a/lib/BoundaryValueDiffEqMIRK/src/interpolation.jl b/lib/BoundaryValueDiffEqMIRK/src/interpolation.jl
index ed806d54..85da8cff 100644
--- a/lib/BoundaryValueDiffEqMIRK/src/interpolation.jl
+++ b/lib/BoundaryValueDiffEqMIRK/src/interpolation.jl
@@ -170,6 +170,18 @@ function (s::EvalSol{C})(tval::Number) where {C <: MIRKCache}
     return z
 end
 
+function (s::EvalSol{C})(tval::Number, T::Type{Val{1}}) where {C <: MIRKCache}
+    (; t, u, cache) = s
+    (; alg, stage, k_discrete, mesh_dt) = cache
+    z′ = zero(last(u))
+    ii = interval(t, tval)
+    dt = mesh_dt[ii]
+    τ = (tval - t[ii]) / dt
+    _, w′ = interp_weights(τ, alg)
+    __maybe_matmul!(z′, k_discrete[ii].du[:, 1:stage], w′[1:stage])
+    return z′
+end
+
 @inline function evalsol_interp_weights(τ::T, ::MIRK2) where {T}
     w = [0, τ * (1 - τ / 2), τ^2 / 2]
 
diff --git a/lib/BoundaryValueDiffEqMIRK/test/mirk_basic_tests.jl b/lib/BoundaryValueDiffEqMIRK/test/mirk_basic_tests.jl
index 70b0c2e4..30249075 100644
--- a/lib/BoundaryValueDiffEqMIRK/test/mirk_basic_tests.jl
+++ b/lib/BoundaryValueDiffEqMIRK/test/mirk_basic_tests.jl
@@ -381,3 +381,50 @@ end
             prob, MIRK4(), dt = 0.05, abstol = 1e-5, controller = error_control)
     end
 end
+
+@testitem "Derivative boundary conditions" begin
+    prob_bvp_nonlinear_14_f(y1, y2, y3, y4) = y2 * y3 - y1 * y4
+    function prob_bvp_nonlinear_14_f!(du, u, p, t)
+        du[1] = u[2]
+        du[2] = u[3]
+        du[3] = u[4]
+        du[4] = prob_bvp_nonlinear_14_f(u[1], u[2], u[3], u[4])
+    end
+    function prob_bvp_nonlinear_14_bc!(res, u, p, t)
+        res[1] = u(0.0)[1]
+        res[2] = u(0.0)[2]
+        res[3] = u(1.0)[1] - 1
+        res[4] = u(1.0)[2]
+    end
+    prob_bvp_nonlinear_14_function = BVPFunction(
+        prob_bvp_nonlinear_14_f!, prob_bvp_nonlinear_14_bc!)
+    prob_bvp_nonlinear_14_tspan = (0.0, 1.0)
+
+    prob_bvp_nonlinear_normal = BVProblem(
+        prob_bvp_nonlinear_14_function, [1.0, 0.0, 0.0, 0.0], prob_bvp_nonlinear_14_tspan)
+
+    @test_nowarn sol1 = solve(
+        prob_bvp_nonlinear_normal, MIRK4(), dt = 0.01, adaptive = false)
+
+    prob_bvp_nonlinear_14_f(y1, y2, y3, y4) = y2 * y3 - y1 * y4
+    function prob_bvp_nonlinear_14_f!(du, u, p, t)
+        du[1] = u[2]
+        du[2] = u[3]
+        du[3] = u[4]
+        du[4] = prob_bvp_nonlinear_14_f(u[1], u[2], u[3], u[4])
+    end
+    function prob_bvp_nonlinear_bc!(res, u, p, t)
+        res[1] = u(0.0)[1]
+        res[2] = u(0.0, Val{1})[1]
+        res[3] = u(1.0)[1] - 1
+        res[4] = u(1.0, Val{1})[1]
+    end
+    prob_bvp_nonlinear_function = BVPFunction(
+        prob_bvp_nonlinear_14_f!, prob_bvp_nonlinear_bc!)
+    prob_bvp_nonlinear_14_tspan = (0.0, 1.0)
+
+    prob_bvp_nonlinear_14 = BVProblem(
+        prob_bvp_nonlinear_function, [1.0, 0.0, 0.0, 0.0], prob_bvp_nonlinear_14_tspan)
+
+    @test_nowarn sol2 = solve(prob_bvp_nonlinear_14, MIRK4(), dt = 0.01, adaptive = false)
+end