diff --git a/lib/NeuralLyapunovProblemLibrary/Project.toml b/lib/NeuralLyapunovProblemLibrary/Project.toml
index 4711cf6..5e96318 100644
--- a/lib/NeuralLyapunovProblemLibrary/Project.toml
+++ b/lib/NeuralLyapunovProblemLibrary/Project.toml
@@ -7,6 +7,7 @@ version = "0.0.1"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 ModelingToolkit = "961ee093-0014-501f-94e3-6117800e7a78"
 Rotations = "6038ab10-8711-5258-84ad-4b1120ba62dc"
+SciMLBase = "0bca4576-84f4-4d90-8ffe-ffa030f20462"
 
 [weakdeps]
 Plots = "91a5bcdd-55d7-5caf-9e0b-520d859cae80"
@@ -22,6 +23,7 @@ OrdinaryDiffEq = "6.91.0"
 Plots = "1.40.9"
 Rotations = "1.7.1"
 SafeTestsets = "0.1.0"
+SciMLBase = "2.75.1"
 StableRNGs = "1.0.2"
 julia = "1.11"
 
diff --git a/lib/NeuralLyapunovProblemLibrary/src/NeuralLyapunovProblemLibrary.jl b/lib/NeuralLyapunovProblemLibrary/src/NeuralLyapunovProblemLibrary.jl
index d3554eb..dd33ac5 100644
--- a/lib/NeuralLyapunovProblemLibrary/src/NeuralLyapunovProblemLibrary.jl
+++ b/lib/NeuralLyapunovProblemLibrary/src/NeuralLyapunovProblemLibrary.jl
@@ -1,6 +1,7 @@
 module NeuralLyapunovProblemLibrary
 
 using ModelingToolkit
+using SciMLBase: NullParameters
 using LinearAlgebra
 using Rotations: RotZXY
 
diff --git a/lib/NeuralLyapunovProblemLibrary/src/double_pendulum.jl b/lib/NeuralLyapunovProblemLibrary/src/double_pendulum.jl
index 1cafe4d..ebe5918 100644
--- a/lib/NeuralLyapunovProblemLibrary/src/double_pendulum.jl
+++ b/lib/NeuralLyapunovProblemLibrary/src/double_pendulum.jl
@@ -1,5 +1,5 @@
 """
-    DoublePendulum(; actuation=:fully_actuated, name)
+    DoublePendulum(; actuation=:fully_actuated, name, defaults)
 
 Create an `ODESystem` representing an undamped double pendulum.
 
@@ -39,8 +39,11 @@ The four actuation modes are described in the table below and selected via `actu
   - `m1`: mass of the first pendulum.
   - `m2`: mass of the second pendulum.
   - `g`: gravitational acceleration (defaults to 9.81).
+
+Users may optionally provide default values of the parameters through `defaults`: a vector
+of the default values for `[I1, I2, l1, l2, lc1, lc2, m1, m2, g]`.
 """
-function DoublePendulum(; actuation=:fully_actuated, name)
+function DoublePendulum(; actuation=:fully_actuated, name, defaults=NullParameters())
     @independent_variables t
     Dt = Differential(t); DDt = Dt^2
 
@@ -62,17 +65,23 @@ function DoublePendulum(; actuation=:fully_actuated, name)
     q = [θ1, θ2]
     params = [I1, I2, l1, l2, lc1, lc2, m1, m2, g]
 
+    kwargs = if defaults == NullParameters()
+        (; name = name)
+    else
+        (; name = name, defaults = Dict(params .=> defaults))
+    end
+
     if actuation == :fully_actuated
         ########################## Fully-actuated double pendulum ##########################
         @variables τ1(t) [input = true] τ2(t) [input = true]
         u = [τ1, τ2]
 
         eqs = DDt.(q) .~ M \ (-C * Dt.(q) + G + u)
-        return ODESystem(eqs, t, vcat(q, u), params; name)
+        return ODESystem(eqs, t, vcat(q, u), params; kwargs...)
     elseif actuation == :undriven
         ############################# Undriven double pendulum #############################
         eqs = DDt.(q) .~ M \ (-C * Dt.(q) + G)
-        return ODESystem(eqs, t, q, params; name)
+        return ODESystem(eqs, t, q, params; kwargs...)
     else
         ########################## Underactuated double pendulum ###########################
         @variables τ(t) [input = true]
@@ -82,13 +91,13 @@ function DoublePendulum(; actuation=:fully_actuated, name)
             B = [0, 1]
             eqs = DDt.(q) .~ M \ (-C * Dt.(q) + G + B * τ)
 
-            return ODESystem(eqs, t, vcat(q, τ), params; name)
+            return ODESystem(eqs, t, vcat(q, τ), params; kwargs...)
         elseif actuation == :pendubot
             ################################### Pendubot ###################################
             B = [1, 0]
             eqs = DDt.(q) .~ M \ (-C * Dt.(q) + G + B * τ)
 
-            return ODESystem(eqs, t, vcat(q, τ), params; name)
+            return ODESystem(eqs, t, vcat(q, τ), params; kwargs...)
         else
             error(
                 "Invalid actuation for DoublePendulum. Received actuation = :",
@@ -99,18 +108,22 @@ function DoublePendulum(; actuation=:fully_actuated, name)
 end
 
 """
-    Acrobot(; name)
+    Acrobot(; name, defaults)
 
-Alias for [`DoublePendulum(; actuation = :acrobot, name)`](@ref).
+Alias for [`DoublePendulum(; actuation = :acrobot, name, defaults)`](@ref).
 """
-Acrobot(; name) = DoublePendulum(; actuation = :acrobot, name)
+function Acrobot(; name, defaults=NullParameters())
+    return DoublePendulum(; actuation = :acrobot, name, defaults)
+end
 
 """
-    Pendubot(; name)
+    Pendubot(; name, defaults)
 
-Alias for [`DoublePendulum(; actuation = :pendubot, name)`](@ref).
+Alias for [`DoublePendulum(; actuation = :pendubot, name, defaults)`](@ref).
 """
-Pendubot(; name) = DoublePendulum(; actuation = :pendubot, name)
+function Pendubot(; name, defaults=NullParameters())
+    DoublePendulum(; actuation = :pendubot, name, defaults)
+end
 
 """
     plot_double_pendulum(θ1, θ2, p, t; title)
diff --git a/lib/NeuralLyapunovProblemLibrary/src/pendulum.jl b/lib/NeuralLyapunovProblemLibrary/src/pendulum.jl
index a38c297..8e81189 100644
--- a/lib/NeuralLyapunovProblemLibrary/src/pendulum.jl
+++ b/lib/NeuralLyapunovProblemLibrary/src/pendulum.jl
@@ -1,5 +1,5 @@
 """
-    Pendulum(; driven = true, name)
+    Pendulum(; driven = true, name, defaults)
 
 Create an `ODESystem` representing a damped, driven or undriven pendulum, depending on the
 value of driven (defaults to `true`, i.e., driven pendulum).
@@ -13,6 +13,9 @@ divided by the moment of inertia of the pendulum around the pivot (`driven = fal
 
 The name of the `ODESystem` is `name`.
 
+Users may optionally provide default values of the parameters through `defaults`: a vector
+of the default values for `[ζ, ω_0]`.
+
 # Example
 
 ```jldoctest; output = false
@@ -35,37 +38,32 @@ sqrt(sum(abs2, [x_end, y_end, ω_end] .- [0, -1, 0])) < 1e-4
 true
 ```
 """
-function Pendulum(; driven = true, name)
+function Pendulum(; driven = true, name, defaults=NullParameters())
     @independent_variables t
     Dt = Differential(t); DDt = Dt^2
 
     @variables θ(t) τ(t) [input = true]
     @parameters ζ ω_0
 
-    if driven
-        ################################# Driven pendulum ##################################
-        eqs = [DDt(θ) + 2ζ * ω_0 * Dt(θ) + ω_0^2 * sin(θ) ~ τ]
-
-        return ODESystem(
-            eqs,
-            t,
-            [θ, τ],
-            [ζ, ω_0];
-            name
-        )
+    params = [ζ, ω_0]
+    kwargs = if defaults == NullParameters()
+        (; name = name)
     else
-        ################################ Undriven pendulum #################################
-        eqs = [DDt(θ) + 2ζ * ω_0 * Dt(θ) + ω_0^2 * sin(θ) ~ 0]
-
-        return ODESystem(
-            eqs,
-            t,
-            [θ],
-            [ζ, ω_0];
-            name
-        )
-        return pendulum_undriven
+        (; name = name, defaults = Dict(params .=> defaults))
     end
+
+    torque = if driven; τ else 0 end
+    variables = if driven; [θ, τ] else [θ] end
+
+    eqs = [DDt(θ) + 2ζ * ω_0 * Dt(θ) + ω_0^2 * sin(θ) ~ torque]
+
+    return ODESystem(
+        eqs,
+        t,
+        variables,
+        params;
+        kwargs...
+    )
 end
 
 """
diff --git a/lib/NeuralLyapunovProblemLibrary/src/quadrotor.jl b/lib/NeuralLyapunovProblemLibrary/src/quadrotor.jl
index fe81533..c93d9b7 100644
--- a/lib/NeuralLyapunovProblemLibrary/src/quadrotor.jl
+++ b/lib/NeuralLyapunovProblemLibrary/src/quadrotor.jl
@@ -1,6 +1,6 @@
 ##################################### Planar quadrotor #####################################
 """
-    QuadrotorPlanar(; name)
+    QuadrotorPlanar(; name, defaults)
 
 Create an `ODESystem` representing a planar approximation of the quadrotor (technically a
 birotor).
@@ -31,8 +31,11 @@ The name of the `ODESystem` is `name`.
   - `g`: gravitational acceleration in the direction of the negative ``y``-axis (defaults to
     9.81).
   - `r`: distance from center of mass to each rotor.
+
+Users may optionally provide default values of the parameters through `defaults`: a vector
+of the default values for `[m, I_quad, g, r]`.
 """
-function QuadrotorPlanar(; name)
+function QuadrotorPlanar(; name, defaults=NullParameters())
     @independent_variables t
     Dt = Differential(t); DDt = Dt^2
     @variables x(t) y(t) θ(t)
@@ -49,13 +52,20 @@ function QuadrotorPlanar(; name)
         I_quad * DDt(θ) ~ r * (ũ1 - ũ2)
     ]
 
-    return ODESystem(eqs, t, [x, y, θ, u1, u2], [m, I_quad, g, r]; name)
+    params = [m, I_quad, g, r]
+    kwargs = if defaults == NullParameters()
+        (; name = name)
+    else
+        (; name = name, defaults = Dict(params .=> defaults))
+    end
+
+    return ODESystem(eqs, t, [x, y, θ, u1, u2], params; kwargs...)
 end
 
 ####################################### 3D quadrotor #######################################
 
 """
-    Quadrotor3D(; name)
+    Quadrotor3D(; name, defaults)
 
 Create an `ODESystem` representing a quadrotor in 3D space.
 
@@ -99,8 +109,11 @@ The model calculates individual rotor thrusts and replaces any negative values w
             I_{xz} & I_{yz} & I_{zz}
         \\end{pmatrix}.
     ```
+
+Users may optionally provide default values of the parameters through `defaults`: a vector
+of the default values for `[m, g, Ixx, Ixy, Ixz, Iyy, Iyz, Izz]`.
 """
-function Quadrotor3D(; name)
+function Quadrotor3D(; name, defaults=NullParameters())
     # Model from "Minimum Snap Trajectory Generation and Control for Quadrotors"
     # https://doi.org/10.1109/ICRA.2011.5980409
     @independent_variables t
@@ -152,12 +165,18 @@ function Quadrotor3D(; name)
                 (τ - angular_velocity_world × (inertia_matrix * angular_velocity_world))
     )
 
+    kwargs = if defaults == NullParameters()
+        (; name = name)
+    else
+        (; name = name, defaults = Dict(params .=> defaults))
+    end
+
     return ODESystem(
         eqs,
         t,
         vcat(position_world, attitude, velocity_world, angular_velocity_world, T, τφ, τθ, τψ),
         params;
-        name
+        kwargs...
     )
 end
 
diff --git a/test/damped_pendulum.jl b/test/damped_pendulum.jl
index 5195968..7c04c3a 100644
--- a/test/damped_pendulum.jl
+++ b/test/damped_pendulum.jl
@@ -26,7 +26,7 @@ eqs = [DDt(θ) + 2ζ * ω_0 * Dt(θ) + ω_0^2 * sin(θ) ~ 0.0]
     t,
     [θ],
     [ζ, ω_0];
-    defaults = defaults
+    defaults
 )
 
 dynamics = structural_simplify(dynamics)