diff --git a/ext/SummationByPartsOperatorsOptimForwardDiffExt.jl b/ext/SummationByPartsOperatorsOptimForwardDiffExt.jl
index 8aa4be77..b4801d41 100644
--- a/ext/SummationByPartsOperatorsOptimForwardDiffExt.jl
+++ b/ext/SummationByPartsOperatorsOptimForwardDiffExt.jl
@@ -12,7 +12,9 @@ function SummationByPartsOperators.function_space_operator(basis_functions, node
                                                            derivative_order = 1, accuracy_order = 0,
                                                            options = Options(g_tol = 1e-14, iterations = 10000)) where {T, SourceOfCoefficients}
 
-    @assert derivative_order == 1 "Only first derivative operators are supported"
+    if derivative_order != 1
+        throw(ArgumentError("Derivative order $derivative_order not implemented."))
+    end
     sort!(nodes)
     weights, D = construct_function_space_operator(basis_functions, nodes, source; options = options)
     return MatrixDerivativeOperator(first(nodes), last(nodes), nodes, weights, D, accuracy_order, source)
diff --git a/test/matrix_operators_test.jl b/test/matrix_operators_test.jl
index 1279de58..f7a31092 100644
--- a/test/matrix_operators_test.jl
+++ b/test/matrix_operators_test.jl
@@ -130,6 +130,8 @@ if VERSION >= v"1.9"
     B[N, N] = 1.0
     let basis_functions = [x -> x^i for i in 0:3]
       D = function_space_operator(basis_functions, nodes, source)
+      # Only first-derivative operators are implemented yet
+      @test_throws ArgumentError function_space_operator(basis_functions, nodes, source; derivative_order = 2)
 
       @test grid(D) ≈ nodes
       @test all(isapprox.(D * ones(N), zeros(N); atol = 1e-13))