diff --git a/docs/make.jl b/docs/make.jl
index ed82c38217..2624bad25d 100644
--- a/docs/make.jl
+++ b/docs/make.jl
@@ -35,8 +35,7 @@ makedocs(
     sitename = "ModelingToolkit.jl",
     authors = "Chris Rackauckas",
     modules = [ModelingToolkitBase, ModelingToolkit, MTKFMIExt],
-    clean = true, doctest = false, linkcheck = true,
-    warnonly = [:docs_block, :missing_docs, :cross_references],
+    clean = true, doctest = true, checkdocs = :exports, linkcheck = true,
     linkcheck_ignore = [
         "https://epubs.siam.org/doi/10.1137/0903023",
         # this link tends to fail linkcheck stochastically and often takes much longer to succeed
diff --git a/docs/src/API/variables.md b/docs/src/API/variables.md
index 291a7ac981..98d5b2b613 100644
--- a/docs/src/API/variables.md
+++ b/docs/src/API/variables.md
@@ -11,6 +11,7 @@ ModelingToolkit to attach additional metadata.
 @parameters
 @constants
 @brownians
+@brownian
 ```
 
 Symbolic variables can have metadata attached to them. The defaults and guesses assigned
diff --git a/docs/src/internals/structural_transformation.md b/docs/src/internals/structural_transformation.md
index 0d4034e475..7da7352852 100644
--- a/docs/src/internals/structural_transformation.md
+++ b/docs/src/internals/structural_transformation.md
@@ -85,4 +85,4 @@ InducedCondensationGraph
 MatchedCondensationGraph
 Unassigned
 unassigned
-```
\ No newline at end of file
+```
diff --git a/lib/ModelingToolkitBase/Project.toml b/lib/ModelingToolkitBase/Project.toml
index b9f79e21ca..58565c0068 100644
--- a/lib/ModelingToolkitBase/Project.toml
+++ b/lib/ModelingToolkitBase/Project.toml
@@ -177,6 +177,7 @@ SafeTestsets = "0.1"
 SciMLBase = "3.44"
 SciMLPublic = "1.0.0"
 SciMLStructures = "1.7"
+SciMLTesting = "2.4"
 Serialization = "1"
 Setfield = "1"
 SimpleNonlinearSolve = "2.11.1"
@@ -231,6 +232,7 @@ REPL = "3fa0cd96-eef1-5676-8a61-b3b8758bbffb"
 Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"
 ReferenceTests = "324d217c-45ce-50fc-942e-d289b448e8cf"
 SafeTestsets = "1bc83da4-3b8d-516f-aca4-4fe02f6d838f"
+SciMLTesting = "09d9d899-5365-40a9-917a-5f67fddea283"
 StableRNGs = "860ef19b-820b-49d6-a774-d7a799459cd3"
 Statistics = "10745b16-79ce-11e8-11f9-7d13ad32a3b2"
 SteadyStateDiffEq = "9672c7b4-1e72-59bd-8a11-6ac3964bc41f"
@@ -239,4 +241,4 @@ Sundials = "c3572dad-4567-51f8-b174-8c6c989267f4"
 Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
 
 [targets]
-test = ["BenchmarkTools", "BoundaryValueDiffEqMIRK", "BoundaryValueDiffEqAscher", "ControlSystemsBase", "DataInterpolations", "DelayDiffEq", "NonlinearSolve", "ForwardDiff", "ModelingToolkitStandardLibrary", "NonlinearSolveBase", "Optimization", "OptimizationIpopt", "OptimizationOptimJL", "OrdinaryDiffEq", "OrdinaryDiffEqDefault", "REPL", "Random", "ReferenceTests", "SafeTestsets", "StableRNGs", "Statistics", "SteadyStateDiffEq", "Test", "StochasticDiffEq", "Sundials", "Pkg", "OrdinaryDiffEqNonlinearSolve", "Logging", "OptimizationBase", "LinearSolve", "Latexify", "Distributed", "DiffEqNoiseProcess", "DynamicQuantities", "OrdinaryDiffEqSDIRK", "OrdinaryDiffEqRosenbrock", "OrdinaryDiffEqBDF", "OrdinaryDiffEqFunctionMap", "OrdinaryDiffEqTsit5"]
+test = ["BenchmarkTools", "BoundaryValueDiffEqMIRK", "BoundaryValueDiffEqAscher", "ControlSystemsBase", "DataInterpolations", "DelayDiffEq", "NonlinearSolve", "ForwardDiff", "ModelingToolkitStandardLibrary", "NonlinearSolveBase", "Optimization", "OptimizationIpopt", "OptimizationOptimJL", "OrdinaryDiffEq", "OrdinaryDiffEqDefault", "REPL", "Random", "ReferenceTests", "SafeTestsets", "SciMLTesting", "StableRNGs", "Statistics", "SteadyStateDiffEq", "Test", "StochasticDiffEq", "Sundials", "Pkg", "OrdinaryDiffEqNonlinearSolve", "Logging", "OptimizationBase", "LinearSolve", "Latexify", "Distributed", "DiffEqNoiseProcess", "DynamicQuantities", "OrdinaryDiffEqSDIRK", "OrdinaryDiffEqRosenbrock", "OrdinaryDiffEqBDF", "OrdinaryDiffEqFunctionMap", "OrdinaryDiffEqTsit5"]
diff --git a/lib/ModelingToolkitBase/src/ModelingToolkitBase.jl b/lib/ModelingToolkitBase/src/ModelingToolkitBase.jl
index 4cad9db4b1..bf0a171214 100644
--- a/lib/ModelingToolkitBase/src/ModelingToolkitBase.jl
+++ b/lib/ModelingToolkitBase/src/ModelingToolkitBase.jl
@@ -155,9 +155,13 @@ abstract type AbstractSystem end
 abstract type IntermediateDeprecationSystem <: AbstractSystem end
 
 """
-    independent_variable(x)
+    independent_variable
 
-Return the independent variable associated with `x`.
+Generic function for querying the primary independent variable of a system-like object.
+
+Most users should call [`independent_variables`](@ref), which returns the independent
+variables as a vector. Packages that define custom system types may extend
+`independent_variable` when a scalar independent-variable interface is required.
 """
 function independent_variable end
 
@@ -255,42 +259,42 @@ include("inputoutput.jl")
 
 include("deprecations.jl")
 
-const t_nounits = let
-    only(@independent_variables t)
-end
-const D_nounits = Differential(t_nounits)
-
-@doc """
-    t
-
-Default independent variable with units.
 """
-t
+    t_nounits
 
-@doc """
-    D
+Unitless default independent variable used by ModelingToolkit examples and constructors.
 
-Default differential operator with respect to [`t`](@ref).
-"""
-D
+# Examples
 
-@doc """
-    t_nounits
+```julia
+using ModelingToolkitBase
 
-Default independent variable without units.
+t = ModelingToolkitBase.t_nounits
+```
 """
-t_nounits
+const t_nounits = let
+    only(@independent_variables t)
+end
 
-@doc """
+"""
     D_nounits
 
-Default differential operator with respect to [`t_nounits`](@ref).
+Default unitless differential operator `Differential(t_nounits)`.
+
+# Examples
+
+```julia
+using ModelingToolkitBase
+
+D = ModelingToolkitBase.D_nounits
+```
 """
-D_nounits
+const D_nounits = Differential(t_nounits)
 
 export CompilerOptions
 export ODEFunction, convert_system_indepvar,
-    System, OptimizationSystem, JumpSystem, SDESystem, NonlinearSystem, ODESystem
+    System, OptimizationSystem, JumpSystem, SDESystem, NonlinearSystem, ODESystem,
+    DiscreteSystem, ImplicitDiscreteSystem
 export SDEFunction
 export DiscreteProblem, DiscreteFunction
 export ImplicitDiscreteProblem, ImplicitDiscreteFunction
diff --git a/lib/ModelingToolkitBase/src/deprecations.jl b/lib/ModelingToolkitBase/src/deprecations.jl
index e880dbebe4..6fc4a9db2f 100644
--- a/lib/ModelingToolkitBase/src/deprecations.jl
+++ b/lib/ModelingToolkitBase/src/deprecations.jl
@@ -1,19 +1,33 @@
-@deprecate structural_simplify(sys; kwargs...) mtkcompile(sys; kwargs...)
-@deprecate structural_simplify(sys, io; kwargs...) mtkcompile(
-    sys; inputs = io[1], outputs = io[2], kwargs...
-)
-
-@doc """
+"""
     structural_simplify(sys; kwargs...)
+    structural_simplify(sys, io; kwargs...)
 
 Deprecated alias for [`mtkcompile`](@ref).
+
+Use `mtkcompile(sys; kwargs...)` in new code.
 """
-structural_simplify
+function structural_simplify(sys; kwargs...)
+    Base.depwarn(
+        "`structural_simplify(sys; kwargs...)` is deprecated, use `mtkcompile(sys; kwargs...)` instead.",
+        :structural_simplify,
+    )
+    return mtkcompile(sys; kwargs...)
+end
+
+function structural_simplify(sys, io; kwargs...)
+    Base.depwarn(
+        "`structural_simplify(sys, io; kwargs...)` is deprecated, use `mtkcompile(sys; kwargs...)` instead.",
+        :structural_simplify,
+    )
+    return mtkcompile(sys; inputs = io[1], outputs = io[2], kwargs...)
+end
 
 """
-    @mtkbuild sys = Constructor(args...; kwargs...)
+    @mtkbuild expr
 
 Deprecated alias for [`@mtkcompile`](@ref).
+
+Use `@mtkcompile` in new code.
 """
 macro mtkbuild(exprs...)
     return quote
@@ -27,14 +41,14 @@ macro mtkbuild(exprs...)
     end |> esc
 end
 
-const ODESystem = IntermediateDeprecationSystem
-
-@doc """
+"""
     ODESystem(args...; kwargs...)
 
-Deprecated alias constructor for [`System`](@ref).
+Deprecated alias for [`System`](@ref).
+
+Use `System(args...; kwargs...)` in new code.
 """
-ODESystem
+const ODESystem = IntermediateDeprecationSystem
 
 function IntermediateDeprecationSystem(args...; kwargs...)
     Base.depwarn(
@@ -45,23 +59,47 @@ function IntermediateDeprecationSystem(args...; kwargs...)
     return System(args...; kwargs...)
 end
 
-for T in [:NonlinearSystem, :DiscreteSystem, :ImplicitDiscreteSystem]
-    @eval @deprecate $T(args...; kwargs...) System(args...; kwargs...)
+"""
+    NonlinearSystem(args...; kwargs...)
+
+Deprecated alias constructor for [`System`](@ref). Use `System(args...; kwargs...)` in
+new code.
+"""
+function NonlinearSystem(args...; kwargs...)
+    Base.depwarn(
+        "`NonlinearSystem(args...; kwargs...)` is deprecated, use `System(args...; kwargs...)` instead.",
+        :NonlinearSystem,
+    )
+    return System(args...; kwargs...)
 end
 
-@doc """
+"""
     DiscreteSystem(args...; kwargs...)
 
-Deprecated alias constructor for [`System`](@ref).
+Deprecated alias constructor for [`System`](@ref). Use `System(args...; kwargs...)` in
+new code.
 """
-DiscreteSystem
+function DiscreteSystem(args...; kwargs...)
+    Base.depwarn(
+        "`DiscreteSystem(args...; kwargs...)` is deprecated, use `System(args...; kwargs...)` instead.",
+        :DiscreteSystem,
+    )
+    return System(args...; kwargs...)
+end
 
-@doc """
+"""
     ImplicitDiscreteSystem(args...; kwargs...)
 
-Deprecated alias constructor for [`System`](@ref).
+Deprecated alias constructor for [`System`](@ref). Use `System(args...; kwargs...)` in
+new code.
 """
-ImplicitDiscreteSystem
+function ImplicitDiscreteSystem(args...; kwargs...)
+    Base.depwarn(
+        "`ImplicitDiscreteSystem(args...; kwargs...)` is deprecated, use `System(args...; kwargs...)` instead.",
+        :ImplicitDiscreteSystem,
+    )
+    return System(args...; kwargs...)
+end
 
 for T in [
         :ODEProblem, :DDEProblem, :SDEProblem, :SDDEProblem, :DAEProblem,
@@ -274,6 +312,8 @@ end
     @brownian xs...
 
 Deprecated alias for [`@brownians`](@ref).
+
+Use `@brownians` in new code.
 """
 macro brownian(xs...)
     return quote
diff --git a/lib/ModelingToolkitBase/src/discretedomain.jl b/lib/ModelingToolkitBase/src/discretedomain.jl
index 4133b3b4d9..959975ff3d 100644
--- a/lib/ModelingToolkitBase/src/discretedomain.jl
+++ b/lib/ModelingToolkitBase/src/discretedomain.jl
@@ -13,10 +13,10 @@ $(FIELDS)
 # Examples
 
 ```jldoctest
-julia> using Symbolics
+julia> using ModelingToolkitBase
 
-julia> Δ = Shift(t)
-(::Shift) (generic function with 2 methods)
+julia> Δ = Shift(ModelingToolkitBase.t_nounits)
+Shift(t, 1)
 ```
 """
 struct Shift <: Operator
@@ -220,12 +220,10 @@ $(FIELDS)
 # Examples
 
 ```jldoctest
-julia> using Symbolics
-
-julia> t = ModelingToolkit.t_nounits
+julia> using ModelingToolkitBase
 
 julia> Δ = Sample(0.01)
-(::Sample) (generic function with 2 methods)
+Sample(SciMLBase.PeriodicClock(0.01, 0.0))
 ```
 """
 struct Sample <: Operator
diff --git a/lib/ModelingToolkitBase/src/problems/bvproblem.jl b/lib/ModelingToolkitBase/src/problems/bvproblem.jl
index 114a66a98f..de6f1cc3e6 100644
--- a/lib/ModelingToolkitBase/src/problems/bvproblem.jl
+++ b/lib/ModelingToolkitBase/src/problems/bvproblem.jl
@@ -13,6 +13,7 @@ function (w::BVPStackedControlRHS)(du, xu, p, t)
 end
 (w::BVPStackedControlRHS)(xu, p, t) = w.f(@view(xu[1:(w.nx)]), @view(xu[(w.nx + 1):end]), p, t)
 
+"""$(problem_docstring(SciMLBase.BVProblem, ODEFunction, true; init = false, extra_body = BV_EXTRA_BODY))"""
 @fallback_iip_specialize function SciMLBase.BVProblem{iip, spec}(
         sys::System, op, tspan;
         check_compatibility = true,
diff --git a/lib/ModelingToolkitBase/src/problems/daeproblem.jl b/lib/ModelingToolkitBase/src/problems/daeproblem.jl
index 8bf07103f4..ce72815365 100644
--- a/lib/ModelingToolkitBase/src/problems/daeproblem.jl
+++ b/lib/ModelingToolkitBase/src/problems/daeproblem.jl
@@ -24,6 +24,7 @@ function generate_DAENLStepData(sys, u0, p, mm, nlstep_compile, nlstep_scc; jac
     )
 end
 
+"""$(function_docstring(DAEFunction, true, [:jac, :tgrad]))"""
 @fallback_iip_specialize function SciMLBase.DAEFunction{iip, spec}(
         sys::System; u0 = nothing, p = nothing, tgrad = false, jac = false,
         t = nothing, eval_expression = false, eval_module = @__MODULE__,
@@ -119,6 +120,7 @@ function SciMLBase.DAEFunction{iip, spec}(
     return maybe_codegen_scimlfn(Val{E}, DAEFunction{iip, spec}, args; kwargs...)
 end
 
+"""$(problem_docstring(SciMLBase.DAEProblem, DAEFunction, true))"""
 @fallback_iip_specialize function SciMLBase.DAEProblem{iip, spec}(
         sys::System, op, tspan;
         callback = nothing, check_length = true, eval_expression = false,
diff --git a/lib/ModelingToolkitBase/src/problems/ddeproblem.jl b/lib/ModelingToolkitBase/src/problems/ddeproblem.jl
index db234fb1bc..1d3d4e8cea 100644
--- a/lib/ModelingToolkitBase/src/problems/ddeproblem.jl
+++ b/lib/ModelingToolkitBase/src/problems/ddeproblem.jl
@@ -1,3 +1,4 @@
+"""$(function_docstring(DDEFunction, true, Symbol[]))"""
 @fallback_iip_specialize function SciMLBase.DDEFunction{iip, spec}(
         sys::System; u0 = nothing, p = nothing, t = nothing, eval_expression = false,
         eval_module = @__MODULE__, expression = Val{false}, checkbounds = false,
@@ -58,6 +59,7 @@ function SciMLBase.DDEFunction{iip, spec}(
     return maybe_codegen_scimlfn(Val{E}, DDEFunction{iip, spec}, args; kwargs...)
 end
 
+"""$(problem_docstring(SciMLBase.DDEProblem, DDEFunction, true))"""
 @fallback_iip_specialize function SciMLBase.DDEProblem{iip, spec}(
         sys::System, op, tspan;
         callback = nothing, check_length = true, checkbounds = false,
diff --git a/lib/ModelingToolkitBase/src/problems/discreteproblem.jl b/lib/ModelingToolkitBase/src/problems/discreteproblem.jl
index c913cd5c88..9a23fdd637 100644
--- a/lib/ModelingToolkitBase/src/problems/discreteproblem.jl
+++ b/lib/ModelingToolkitBase/src/problems/discreteproblem.jl
@@ -1,3 +1,4 @@
+"""$(function_docstring(DiscreteFunction, true, Symbol[]))"""
 @fallback_iip_specialize function SciMLBase.DiscreteFunction{iip, spec}(
         sys::System; u0 = nothing, p = nothing, t = nothing,
         eval_expression = false, eval_module = @__MODULE__, expression = Val{false},
@@ -54,6 +55,7 @@ function SciMLBase.DiscreteFunction{iip, spec}(
     return maybe_codegen_scimlfn(Val{E}, DiscreteFunction{iip, spec}, args; kwargs...)
 end
 
+"""$(problem_docstring(SciMLBase.DiscreteProblem, DiscreteFunction, true))"""
 @fallback_iip_specialize function SciMLBase.DiscreteProblem{iip, spec}(
         sys::System, op, tspan;
         check_compatibility = true, expression = Val{false}, kwargs...
diff --git a/lib/ModelingToolkitBase/src/problems/docs.jl b/lib/ModelingToolkitBase/src/problems/docs.jl
index 937b63ef67..e1266284c5 100644
--- a/lib/ModelingToolkitBase/src/problems/docs.jl
+++ b/lib/ModelingToolkitBase/src/problems/docs.jl
@@ -195,31 +195,6 @@ If the `System` has algebraic equations, like `x(t)^2 + y(t)^2`, the resulting
 `BVProblem` must be solved using BVDAE solvers, such as Ascher.
 """
 
-for (mod, prob, func, istd, kws) in [
-        (SciMLBase, :ODEProblem, ODEFunction, true, (;)),
-        (SciMLBase, :SteadyStateProblem, ODEFunction, false, (;)),
-        (
-            SciMLBase, :BVProblem, ODEFunction, true,
-            (; init = false, extra_body = BV_EXTRA_BODY),
-        ),
-        (SciMLBase, :DAEProblem, DAEFunction, true, (;)),
-        (SciMLBase, :DDEProblem, DDEFunction, true, (;)),
-        (SciMLBase, :SDEProblem, SDEFunction, true, (;)),
-        (SciMLBase, :SDDEProblem, SDDEFunction, true, (;)),
-        (JumpProcesses, :JumpProblem, "inner SciMLFunction", true, (; init = false)),
-        (SciMLBase, :DiscreteProblem, DiscreteFunction, true, (;)),
-        (SciMLBase, :ImplicitDiscreteProblem, ImplicitDiscreteFunction, true, (;)),
-        (SciMLBase, :NonlinearProblem, NonlinearFunction, false, (;)),
-        (SciMLBase, :NonlinearLeastSquaresProblem, NonlinearFunction, false, (;)),
-        (SciMLBase, :OptimizationProblem, OptimizationFunction, false, (; init = false)),
-    ]
-    kwexpr = Expr(:parameters)
-    for (k, v) in pairs(kws)
-        push!(kwexpr.args, Expr(:kw, k, v))
-    end
-    @eval @doc problem_docstring($kwexpr, $mod.$prob, $func, $istd) $mod.$prob
-end
-
 function function_docstring(
         func, istd, optionals; extra_body = "", extra_kwargs = "", extra_kwargs_desc = ""
     )
@@ -347,103 +322,3 @@ function process_optional_function_kwargs(choices::Vector{Symbol})
     end
     return join(map(Base.Fix1(getindex, OPTIONAL_FN_KWARGS_DICT), choices), "\n")
 end
-
-for (mod, func, istd, optionals, kws) in [
-        (SciMLBase, :ODEFunction, true, [:jac, :tgrad], (;)),
-        (SciMLBase, :ODEInputFunction, true, [:inputfn, :jac, :tgrad, :controljac], (;)),
-        (SciMLBase, :DAEFunction, true, [:jac, :tgrad], (;)),
-        (SciMLBase, :DDEFunction, true, Symbol[], (;)),
-        (SciMLBase, :SDEFunction, true, [:jac, :tgrad], (;)),
-        (SciMLBase, :SDDEFunction, true, Symbol[], (;)),
-        (SciMLBase, :DiscreteFunction, true, Symbol[], (;)),
-        (SciMLBase, :ImplicitDiscreteFunction, true, Symbol[], (;)),
-        (SciMLBase, :NonlinearFunction, false, [:resid_prototype, :jac], (;)),
-        (SciMLBase, :IntervalNonlinearFunction, false, Symbol[], (;)),
-        (SciMLBase, :OptimizationFunction, false, [:jac, :grad, :hess, :cons_h, :cons_j], (;)),
-    ]
-    kwexpr = Expr(:parameters)
-    for (k, v) in pairs(kws)
-        push!(kwexpr.args, Expr(:kw, k, v))
-    end
-    @eval @doc function_docstring($kwexpr, $mod.$func, $istd, $optionals) $mod.$func
-end
-
-@doc """
-    SciMLBase.HomotopyNonlinearFunction(sys::System; kwargs...)
-    SciMLBase.HomotopyNonlinearFunction{iip}(sys::System; kwargs...)
-    SciMLBase.HomotopyNonlinearFunction{iip, specialize}(sys::System; kwargs...)
-
-Create a `HomotopyNonlinearFunction` from the given `sys`. `iip` is a boolean indicating
-whether the function should be in-place. `specialization` is a `SciMLBase.AbstractSpecalize`
-subtype indicating the level of specialization of the $func.
-
-# Keyword arguments
-
-- `u0`: The `u0` vector for the corresponding problem, if available. Can be obtained
-  using [`ModelingToolkitBase.get_u0`](@ref).
-- `p`: The parameter object for the corresponding problem, if available. Can be obtained
-  using [`ModelingToolkitBase.get_p`](@ref).
-$EVAL_EXPR_MOD_KWARGS
-- `checkbounds`: Whether to enable bounds checking in the generated code.
-- `simplify`: Whether to `simplify` any symbolically computed jacobians/hessians/etc.
-  This typically improves performance of the generated code but reduces readability.
-- `fraction_cancel_fn`: The function to use to simplify fractions in the polynomial
-  expression. A more powerful function can increase processing time but be able to
-  eliminate more rational functions, thus improving solve time. Should be a function that
-  takes a symbolic expression containing zero or more fraction expressions and returns the
-  simplified expression. While this defaults to `SymbolicUtils.simplify_fractions`, a viable
-  alternative is `SymbolicUtils.quick_cancel`
-
-All keyword arguments are forwarded to the wrapped `NonlinearFunction` constructor.
-""" SciMLBase.HomotopyNonlinearFunction
-
-@doc """
-    SciMLBase.IntervalNonlinearProblem(sys::System, uspan::NTuple{2}, parammap = SciMLBase.NullParameters(); kwargs...)
-
-Create an `IntervalNonlinearProblem` from the given `sys`. This is only valid for a system
-of nonlinear equations with a single equation and unknown. `uspan` is the interval in which
-the root is to be found, and `parammap` is an iterable collection of key-value pairs
-providing values for the parameters in the system.
-
-$TIME_INDEPENDENT_INIT
-
-# Keyword arguments
-
-$PROBLEM_KWARGS
-$(prob_fun_common_kwargs(IntervalNonlinearProblem, false))
-
-All other keyword arguments are forwarded to the `IntervalNonlinearFunction` constructor.
-
-$PROBLEM_INTERNALS_HEADER
-
-$PROBLEM_INTERNAL_KWARGS
-""" SciMLBase.IntervalNonlinearProblem
-
-@doc """
-    SciMLBase.LinearProblem(sys::System, op; kwargs...)
-    SciMLBase.LinearProblem{iip}(sys::System, op; kwargs...)
-
-Build a `LinearProblem` given a system `sys` and operating point `op`. `iip` is a boolean
-indicating whether the problem should be in-place. The operating point should be an
-iterable collection of key-value pairs mapping variables/parameters in the system to the
-(initial) values they should take in `LinearProblem`. Any values not provided will
-fallback to the corresponding default (if present).
-
-Note that since `u0` is optional for `LinearProblem`, values of unknowns do not need to be
-specified in `op` to create a `LinearProblem`. In such a case, `prob.u0` will be `nothing`
-and attempting to symbolically index the problem with an unknown, observable, or expression
-depending on unknowns/observables will error.
-
-Updating the parameters automatically updates the `A` and `b` arrays.
-
-# Keyword arguments
-
-$PROBLEM_KWARGS
-$(prob_fun_common_kwargs(LinearProblem, false))
-
-All other keyword arguments are forwarded to the $func constructor.
-
-$PROBLEM_INTERNALS_HEADER
-
-$PROBLEM_INTERNAL_KWARGS
-""" SciMLBase.LinearProblem
diff --git a/lib/ModelingToolkitBase/src/problems/implicitdiscreteproblem.jl b/lib/ModelingToolkitBase/src/problems/implicitdiscreteproblem.jl
index b941ea3ba5..56495813c4 100644
--- a/lib/ModelingToolkitBase/src/problems/implicitdiscreteproblem.jl
+++ b/lib/ModelingToolkitBase/src/problems/implicitdiscreteproblem.jl
@@ -1,3 +1,4 @@
+"""$(function_docstring(ImplicitDiscreteFunction, true, Symbol[]))"""
 @fallback_iip_specialize function SciMLBase.ImplicitDiscreteFunction{iip, spec}(
         sys::System; u0 = nothing, p = nothing, t = nothing, eval_expression = false,
         eval_module = @__MODULE__, expression = Val{false},
@@ -63,6 +64,7 @@ function SciMLBase.ImplicitDiscreteFunction{iip, spec}(
     )
 end
 
+"""$(problem_docstring(SciMLBase.ImplicitDiscreteProblem, ImplicitDiscreteFunction, true))"""
 @fallback_iip_specialize function SciMLBase.ImplicitDiscreteProblem{iip, spec}(
         sys::System, op, tspan;
         check_compatibility = true, expression = Val{false}, kwargs...
diff --git a/lib/ModelingToolkitBase/src/problems/initializationproblem.jl b/lib/ModelingToolkitBase/src/problems/initializationproblem.jl
index 9418fed0b8..35d1aa2439 100644
--- a/lib/ModelingToolkitBase/src/problems/initializationproblem.jl
+++ b/lib/ModelingToolkitBase/src/problems/initializationproblem.jl
@@ -1,6 +1,4 @@
-struct InitializationProblem{iip, specialization} end
-
-@doc """
+"""
     InitializationProblem(sys::AbstractSystem, t, op = Dict(); kwargs...)
     InitializationProblem{iip}(sys::AbstractSystem, t, op = Dict(); kwargs...)
     InitializationProblem{iip, specialize}(sys::AbstractSystem, t, op = Dict(); kwargs...)
@@ -20,7 +18,8 @@ $INITIALIZEPROB_KWARGS
 $INTERNAL_INITIALIZEPROB_KWARGS
 
 All other keyword arguments are forwarded to the wrapped problem constructor.
-""" InitializationProblem
+"""
+struct InitializationProblem{iip, specialization} end
 
 """
     $(TYPEDSIGNATURES)
diff --git a/lib/ModelingToolkitBase/src/problems/intervalnonlinearproblem.jl b/lib/ModelingToolkitBase/src/problems/intervalnonlinearproblem.jl
index 94965a7cb5..a3dcf27df7 100644
--- a/lib/ModelingToolkitBase/src/problems/intervalnonlinearproblem.jl
+++ b/lib/ModelingToolkitBase/src/problems/intervalnonlinearproblem.jl
@@ -1,3 +1,4 @@
+"""$(function_docstring(IntervalNonlinearFunction, false, Symbol[]))"""
 function SciMLBase.IntervalNonlinearFunction(
         sys::System; u0 = nothing, p = nothing, t = nothing, eval_expression = false,
         eval_module = @__MODULE__, expression = Val{false}, checkbounds = false,
@@ -43,6 +44,7 @@ function SciMLBase.IntervalNonlinearFunction(sys::System, opts::SciMLFunctionOpt
     )
 end
 
+"""$(problem_docstring(SciMLBase.IntervalNonlinearProblem, IntervalNonlinearFunction, false))"""
 function SciMLBase.IntervalNonlinearProblem(
         sys::System, uspan::NTuple{2}, parammap = SciMLBase.NullParameters();
         check_compatibility = true, expression = Val{false}, kwargs...
diff --git a/lib/ModelingToolkitBase/src/problems/jumpproblem.jl b/lib/ModelingToolkitBase/src/problems/jumpproblem.jl
index df8fb04a5e..6b9b0a9be6 100644
--- a/lib/ModelingToolkitBase/src/problems/jumpproblem.jl
+++ b/lib/ModelingToolkitBase/src/problems/jumpproblem.jl
@@ -1,3 +1,4 @@
+"""$(problem_docstring(JumpProcesses.JumpProblem, "inner SciMLFunction", true; init = false))"""
 @fallback_iip_specialize function JumpProcesses.JumpProblem{iip, spec}(
         sys::System, op, tspan::Union{Tuple, Nothing};
         check_compatibility = true, eval_expression = false, eval_module = @__MODULE__,
diff --git a/lib/ModelingToolkitBase/src/problems/linearproblem.jl b/lib/ModelingToolkitBase/src/problems/linearproblem.jl
index 529030964e..c726e72e1a 100644
--- a/lib/ModelingToolkitBase/src/problems/linearproblem.jl
+++ b/lib/ModelingToolkitBase/src/problems/linearproblem.jl
@@ -84,6 +84,7 @@ function LinearFunction{iip}(
     return LinearFunction{iip}(symbolic_interface, A, b)
 end
 
+"""$(problem_docstring(SciMLBase.LinearProblem, LinearFunction, false))"""
 function SciMLBase.LinearProblem(sys::System, op; kwargs...)
     return SciMLBase.LinearProblem{true}(sys, op; kwargs...)
 end
diff --git a/lib/ModelingToolkitBase/src/problems/nonlinearproblem.jl b/lib/ModelingToolkitBase/src/problems/nonlinearproblem.jl
index 715e0668d9..631eebb2bf 100644
--- a/lib/ModelingToolkitBase/src/problems/nonlinearproblem.jl
+++ b/lib/ModelingToolkitBase/src/problems/nonlinearproblem.jl
@@ -1,3 +1,4 @@
+"""$(function_docstring(NonlinearFunction, false, [:resid_prototype, :jac]))"""
 @fallback_iip_specialize function SciMLBase.NonlinearFunction{iip, spec}(
         sys::System; u0 = nothing, p = nothing, t = nothing, jac = false,
         eval_expression = false, eval_module = @__MODULE__, sparse = false,
@@ -113,6 +114,7 @@ function generate_nonlinear_bounds(sys::AbstractSystem, op)
     return lb, ub
 end
 
+"""$(problem_docstring(SciMLBase.NonlinearProblem, NonlinearFunction, false))"""
 @fallback_iip_specialize function SciMLBase.NonlinearProblem{iip, spec}(
         sys::System, op; expression = Val{false}, lb = nothing, ub = nothing,
         check_length = true, check_compatibility = true, kwargs...
@@ -168,6 +170,7 @@ function SciMLBase.AbstractNonlinearProblem(sys::System, op; kwargs...)
     return get_nonlinear_problem_type(sys)(sys, op; kwargs...)
 end
 
+"""$(problem_docstring(SciMLBase.NonlinearLeastSquaresProblem, NonlinearFunction, false))"""
 @fallback_iip_specialize function SciMLBase.NonlinearLeastSquaresProblem{iip, spec}(
         sys::System, op; check_length = false, lb = nothing, ub = nothing,
         check_compatibility = true, expression = Val{false}, kwargs...
diff --git a/lib/ModelingToolkitBase/src/problems/odeproblem.jl b/lib/ModelingToolkitBase/src/problems/odeproblem.jl
index 887f0f8367..171ba0d998 100644
--- a/lib/ModelingToolkitBase/src/problems/odeproblem.jl
+++ b/lib/ModelingToolkitBase/src/problems/odeproblem.jl
@@ -18,6 +18,7 @@ function generate_ODENLStepData(sys, u0, p, mm, nlstep_compile, nlstep_scc; jac
     )
 end
 
+"""$(function_docstring(ODEFunction, true, [:jac, :tgrad]))"""
 Base.@nospecializeinfer @fallback_iip_specialize function SciMLBase.ODEFunction{iip, spec}(
         sys::System; @nospecialize(u0 = nothing), @nospecialize(p = nothing), t = nothing,
         tgrad = false, jac = false,
@@ -149,6 +150,7 @@ Base.@nospecializeinfer function _ode_problem(
     return maybe_codegen_scimlproblem(expression, ODEProblem{_iip}, args; kwargs...)
 end
 
+"""$(problem_docstring(SciMLBase.ODEProblem, ODEFunction, true))"""
 Base.@nospecializeinfer @fallback_iip_specialize function SciMLBase.ODEProblem{iip, spec}(
         sys::System, @nospecialize(op), tspan;
         @nospecialize(callback = nothing), check_length = true, eval_expression = false,
@@ -178,6 +180,7 @@ Base.@nospecializeinfer function SciMLBase.ODEProblem{
     )
 end
 
+"""$(problem_docstring(DiffEqBase.SteadyStateProblem, ODEFunction, false))"""
 @fallback_iip_specialize function DiffEqBase.SteadyStateProblem{iip, spec}(
         sys::System, op; check_length = true, check_compatibility = true,
         expression = Val{false}, kwargs...
diff --git a/lib/ModelingToolkitBase/src/problems/optimizationproblem.jl b/lib/ModelingToolkitBase/src/problems/optimizationproblem.jl
index 3666dc923f..1b87bb9218 100644
--- a/lib/ModelingToolkitBase/src/problems/optimizationproblem.jl
+++ b/lib/ModelingToolkitBase/src/problems/optimizationproblem.jl
@@ -1,3 +1,4 @@
+"""$(function_docstring(OptimizationFunction, false, [:jac, :grad, :hess, :cons_h, :cons_j]))"""
 function SciMLBase.OptimizationFunction(sys::System, args...; kwargs...)
     return OptimizationFunction{true}(sys, args...; kwargs...)
 end
@@ -104,6 +105,7 @@ function SciMLBase.OptimizationFunction{iip}(
     return maybe_codegen_scimlfn(Val{E}, OptimizationFunction{iip}, args; kwargs...)
 end
 
+"""$(problem_docstring(SciMLBase.OptimizationProblem, OptimizationFunction, false; init = false))"""
 function SciMLBase.OptimizationProblem(sys::System, args...; kwargs...)
     return OptimizationProblem{true}(sys, args...; kwargs...)
 end
diff --git a/lib/ModelingToolkitBase/src/problems/sddeproblem.jl b/lib/ModelingToolkitBase/src/problems/sddeproblem.jl
index 1fbe1be5a4..3543714ded 100644
--- a/lib/ModelingToolkitBase/src/problems/sddeproblem.jl
+++ b/lib/ModelingToolkitBase/src/problems/sddeproblem.jl
@@ -1,3 +1,4 @@
+"""$(function_docstring(SDDEFunction, true, Symbol[]))"""
 @fallback_iip_specialize function SciMLBase.SDDEFunction{iip, spec}(
         sys::System; u0 = nothing, p = nothing, t = nothing, expression = Val{false},
         eval_expression = false, eval_module = @__MODULE__,
@@ -60,6 +61,7 @@ function SciMLBase.SDDEFunction{iip, spec}(
     return maybe_codegen_scimlfn(Val{E}, SDDEFunction{iip, spec}, args; kwargs...)
 end
 
+"""$(problem_docstring(SciMLBase.SDDEProblem, SDDEFunction, true))"""
 @fallback_iip_specialize function SciMLBase.SDDEProblem{iip, spec}(
         sys::System, op, tspan;
         callback = nothing, check_length = true, checkbounds = false,
diff --git a/lib/ModelingToolkitBase/src/problems/sdeproblem.jl b/lib/ModelingToolkitBase/src/problems/sdeproblem.jl
index 98cf7a5d3a..8dd376cc21 100644
--- a/lib/ModelingToolkitBase/src/problems/sdeproblem.jl
+++ b/lib/ModelingToolkitBase/src/problems/sdeproblem.jl
@@ -1,3 +1,4 @@
+"""$(function_docstring(SDEFunction, true, [:jac, :tgrad]))"""
 @fallback_iip_specialize function SciMLBase.SDEFunction{iip, spec}(
         sys::System; u0 = nothing, p = nothing, tgrad = false, jac = false,
         t = nothing, eval_expression = false, eval_module = @__MODULE__,
@@ -82,6 +83,7 @@ function SciMLBase.SDEFunction{iip, spec}(
     return maybe_codegen_scimlfn(Val{E}, SDEFunction{iip, spec}, args; kwargs...)
 end
 
+"""$(problem_docstring(SciMLBase.SDEProblem, SDEFunction, true))"""
 @fallback_iip_specialize function SciMLBase.SDEProblem{iip, spec}(
         sys::System, op, tspan;
         callback = nothing, check_length = true, eval_expression = false,
diff --git a/lib/ModelingToolkitBase/src/systems/abstractsystem.jl b/lib/ModelingToolkitBase/src/systems/abstractsystem.jl
index 4522152bf9..bf09e50186 100644
--- a/lib/ModelingToolkitBase/src/systems/abstractsystem.jl
+++ b/lib/ModelingToolkitBase/src/systems/abstractsystem.jl
@@ -1893,9 +1893,32 @@ function state_priorities(sys::AbstractSystem)
 end
 
 """
-    $TYPEDSIGNATURES
+    irreducibles(sys::AbstractSystem)
+
+Return the variables in `sys` and its subsystems that are marked as irreducible.
+
+Irreducible variables are preserved as unknowns during simplification instead of being
+eliminated as observed variables when possible.
+
+# Arguments
+
+- `sys`: system to inspect recursively.
+
+# Returns
+
+An atomic set of symbolic variables, with subsystem variables namespaced into `sys`.
+
+# Examples
+
+```julia
+using ModelingToolkitBase
+using ModelingToolkitBase: t_nounits as t, D_nounits as D
+
+@variables x(t) [irreducible = true]
+@named sys = System([D(x) ~ -x], t)
 
-Get the irreducible variables of a system `sys` and its subsystems.
+irreducibles(sys)
+```
 """
 function irreducibles(sys::AbstractSystem)
     ircs = get_irreducibles(sys)
@@ -1909,9 +1932,20 @@ function irreducibles(sys::AbstractSystem)
 end
 
 """
-    $TYPEDSIGNATURES
+    maybe_zeros(sys::AbstractSystem)
+
+Return variables in `sys` and its subsystems that simplification may constrain to zero.
+
+This is primarily used by structural simplification to track variables introduced or
+retained while handling alias equations.
+
+# Arguments
+
+- `sys`: system to inspect recursively.
+
+# Returns
 
-Get the variables that should be treated as possible zeros in `sys` and its subsystems.
+An atomic set of symbolic variables, with subsystem variables namespaced into `sys`.
 """
 function maybe_zeros(sys::AbstractSystem)
     dds = get_maybe_zeros(sys)
diff --git a/lib/ModelingToolkitBase/src/systems/nonlinear/homotopy_continuation.jl b/lib/ModelingToolkitBase/src/systems/nonlinear/homotopy_continuation.jl
index f02dd95cc5..8c908f2ddc 100644
--- a/lib/ModelingToolkitBase/src/systems/nonlinear/homotopy_continuation.jl
+++ b/lib/ModelingToolkitBase/src/systems/nonlinear/homotopy_continuation.jl
@@ -490,6 +490,7 @@ function handle_rational_polynomials(x, wrt; fraction_cancel_fn = simplify_fract
     return num, den
 end
 
+"""$(function_docstring(HomotopyNonlinearFunction, false, Symbol[]; extra_kwargs = "- `fraction_cancel_fn`: Function used to simplify fractions in polynomial expressions."))"""
 @fallback_iip_specialize function SciMLBase.HomotopyNonlinearFunction{iip, specialize}(
         sys::System; eval_expression = false, eval_module = @__MODULE__,
         p = nothing, fraction_cancel_fn = SymbolicUtils.simplify_fractions,
@@ -539,11 +540,23 @@ end
     )
 end
 
-struct HomotopyContinuationProblem{iip, specialization} end
+"""
+    HomotopyContinuationProblem(sys::System, args...; kwargs...)
+
+Construct a homotopy-continuation problem from a nonlinear `sys`. This problem type is
+intended for solvers that track solution paths from a start system to the target system.
+
+# Example
 
-@doc problem_docstring(
-    HomotopyContinuationProblem, HomotopyNonlinearFunction, false; init = false
-) HomotopyContinuationProblem
+```julia
+@variables x = 1.0
+@mtkcompile sys = System([x^2 - 1 ~ 0])
+prob = HomotopyContinuationProblem(sys, [])
+```
+
+See [`NonlinearProblem`](@ref) for the standard nonlinear-problem constructor.
+"""
+struct HomotopyContinuationProblem{iip, specialization} end
 
 function HomotopyContinuationProblem(sys::System, args...; kwargs...)
     return HomotopyContinuationProblem{true}(sys, args...; kwargs...)
diff --git a/lib/ModelingToolkitBase/src/systems/optimal_control_interface.jl b/lib/ModelingToolkitBase/src/systems/optimal_control_interface.jl
index 3ae621505d..941e023b33 100644
--- a/lib/ModelingToolkitBase/src/systems/optimal_control_interface.jl
+++ b/lib/ModelingToolkitBase/src/systems/optimal_control_interface.jl
@@ -1,15 +1,24 @@
 """
     AbstractCollocation
 
-Abstract supertype for dynamic optimization collocation solver descriptors.
+Abstract supertype for dynamic-optimization collocation method descriptors.
+
+Concrete collocation types are provided by backend extensions such as JuMP, InfiniteOpt,
+CasADi, and Pyomo.
 """
 abstract type AbstractCollocation end
 
 """
-    DynamicOptSolution
+    DynamicOptSolution(model, sol, input_sol)
+
+Solution wrapper returned by dynamic-optimization backends.
+
+# Fields
 
-Container returned by dynamic optimization solves, holding the optimized model, state
-trajectory solution, and optional input trajectory solution.
+- `model`: backend-specific optimization model.
+- `sol`: state trajectory as an `ODESolution`.
+- `input_sol`: controller/input trajectory as an `ODESolution`, or `nothing` when the
+  backend does not return one.
 """
 struct DynamicOptSolution
     model::Any
@@ -143,6 +152,7 @@ end
 
 is_explicit(tableau) = tableau isa DiffEqBase.ExplicitRKTableau
 
+"""$(function_docstring(ODEInputFunction, true, [:inputfn, :jac, :tgrad, :controljac]))"""
 @fallback_iip_specialize function SciMLBase.ODEInputFunction{iip, specialize}(
         sys::System;
         inputs = default_codegen_inputs(sys),
diff --git a/lib/ModelingToolkitBase/src/systems/unit_check.jl b/lib/ModelingToolkitBase/src/systems/unit_check.jl
index 8e82bf0d09..738c27c3ca 100644
--- a/lib/ModelingToolkitBase/src/systems/unit_check.jl
+++ b/lib/ModelingToolkitBase/src/systems/unit_check.jl
@@ -11,6 +11,14 @@ end
 struct PleaseImportDynamicQuantities end
 global t::Union{PleaseImportDynamicQuantities, Num} = PleaseImportDynamicQuantities()
 
+@doc """
+    t
+
+Default independent variable with units. Only usable once DynamicQuantities.jl is
+loaded; see [`t_nounits`](@ref) for the unitless default.
+"""
+t
+
 function Base.show(io::IO, ::PleaseImportDynamicQuantities)
     return __import_dynamic_quantities()
 end
@@ -23,3 +31,11 @@ function __import_dynamic_quantities(_...)
     )
 end
 global D::Union{typeof(__import_dynamic_quantities), Differential} = __import_dynamic_quantities
+
+@doc """
+    D
+
+Default differential operator with respect to [`t`](@ref). Only usable once
+DynamicQuantities.jl is loaded; see [`D_nounits`](@ref) for the unitless default.
+"""
+D
diff --git a/lib/ModelingToolkitBase/src/utils.jl b/lib/ModelingToolkitBase/src/utils.jl
index cab12a4d94..b25bc90634 100644
--- a/lib/ModelingToolkitBase/src/utils.jl
+++ b/lib/ModelingToolkitBase/src/utils.jl
@@ -580,9 +580,18 @@ function collect_guesses!(guesses::SymmapT, vars::Vector{SymbolicT})
 end
 
 """
-    $TYPEDSIGNATURES
+    collect_var_to_name!(vars::Dict{Symbol, SymbolicT}, xs::Vector{SymbolicT})
+
+Populate `vars` with mappings from symbolic variable names to variables.
+
+# Arguments
+
+- `vars`: dictionary updated in place.
+- `xs`: symbolic variables to inspect.
+
+# Returns
 
-Populate `vars` with a mapping from the name of each symbolic variable in `xs` to that variable.
+`nothing`. Throws `ArgumentError` if two distinct variables have the same name.
 """
 function collect_var_to_name!(vars::Dict{Symbol, SymbolicT}, xs::Vector{SymbolicT})
     for x in xs
diff --git a/lib/ModelingToolkitBase/src/variables.jl b/lib/ModelingToolkitBase/src/variables.jl
index 7493deb163..2574e4a856 100644
--- a/lib/ModelingToolkitBase/src/variables.jl
+++ b/lib/ModelingToolkitBase/src/variables.jl
@@ -504,7 +504,28 @@ end
 """
     setnominal(x, val)
 
-Attach nominal value `val` to symbolic variable `x`.
+Return `x` with nominal-value metadata set to `val`.
+
+# Arguments
+
+- `x`: symbolic variable to annotate.
+- `val`: nominal value used for scaling and numerical conditioning.
+
+# Returns
+
+A symbolic variable equivalent to `x` with updated `VariableNominal` metadata.
+
+# Examples
+
+```julia
+using ModelingToolkitBase
+
+@variables x
+x = setnominal(x, 10.0)
+getnominal(x)
+```
+
+See also [`getnominal`](@ref) and [`hasnominal`](@ref).
 """
 function setnominal(x::Num, val)
     return setmetadata(x, VariableNominal, val)
@@ -716,6 +737,14 @@ isbrownian(s) = getvariabletype(s) === BROWNIAN
 $(SIGNATURES)
 
 Define one or more Brownian variables.
+
+# Examples
+
+```julia
+using ModelingToolkitBase
+
+@brownians B
+```
 """
 macro brownians(xs...)
     all(
diff --git a/lib/ModelingToolkitBase/test/extensions/Project.toml b/lib/ModelingToolkitBase/test/extensions/Project.toml
index 7b9d389b13..e3f3fa50e3 100644
--- a/lib/ModelingToolkitBase/test/extensions/Project.toml
+++ b/lib/ModelingToolkitBase/test/extensions/Project.toml
@@ -27,5 +27,5 @@ ModelingToolkitBase = {path = "../.."}
 
 [compat]
 SafeTestsets = "0.1, 1"
-SciMLTesting = "1, 2.1"
+SciMLTesting = "2.8"
 TaylorDiff = "0.3.5"
diff --git a/lib/ModelingToolkitBase/test/optimization/Project.toml b/lib/ModelingToolkitBase/test/optimization/Project.toml
index 525728e9f6..cfb9b2aef8 100644
--- a/lib/ModelingToolkitBase/test/optimization/Project.toml
+++ b/lib/ModelingToolkitBase/test/optimization/Project.toml
@@ -38,4 +38,4 @@ OrdinaryDiffEqExplicitTableaus = "2"
 OrdinaryDiffEqImplicitTableaus = "2"
 Pyomo = "2.0.2"
 SafeTestsets = "0.1, 1"
-SciMLTesting = "1, 2.1"
+SciMLTesting = "2.8"
diff --git a/lib/ModelingToolkitBase/test/qa/Project.toml b/lib/ModelingToolkitBase/test/qa/Project.toml
index e5ef2d9282..ef11d67946 100644
--- a/lib/ModelingToolkitBase/test/qa/Project.toml
+++ b/lib/ModelingToolkitBase/test/qa/Project.toml
@@ -4,7 +4,11 @@ JET = "c3a54625-cd67-489e-a8e7-0a5a0ff4e31b"
 ModelingToolkitBase = "7771a370-6774-4173-bd38-47e70ca0b839"
 SciMLBase = "0bca4576-84f4-4d90-8ffe-ffa030f20462"
 SciMLStructures = "53ae85a6-f571-4167-b2af-e1d143709226"
+SciMLTesting = "09d9d899-5365-40a9-917a-5f67fddea283"
 Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
 
 [sources]
 ModelingToolkitBase = {path = "../.."}
+
+[compat]
+SciMLTesting = "2.4"
diff --git a/lib/ModelingToolkitBase/test/qa/aqua.jl b/lib/ModelingToolkitBase/test/qa/aqua.jl
index 7bbb4d629d..9fa5ece7d8 100644
--- a/lib/ModelingToolkitBase/test/qa/aqua.jl
+++ b/lib/ModelingToolkitBase/test/qa/aqua.jl
@@ -1,6 +1,6 @@
 using ModelingToolkitBase
 using Aqua
+using JET
+using SciMLTesting
 
-# This is loaded only by MTKBifurcationKitExt, while remaining a hard dependency so
-# loading BifurcationKit alone still activates the extension.
-Aqua.test_all(ModelingToolkitBase; stale_deps = (; ignore = [:SimpleNonlinearSolve]))
+run_qa(ModelingToolkitBase; Aqua, JET, jet = true)
diff --git a/src/ModelingToolkit.jl b/src/ModelingToolkit.jl
index 896b317e6f..c372bc511e 100644
--- a/src/ModelingToolkit.jl
+++ b/src/ModelingToolkit.jl
@@ -175,6 +175,38 @@ export get_sensitivity_function, get_comp_sensitivity_function,
     get_looptransfer_function, get_sensitivity, get_comp_sensitivity, get_looptransfer
 export isolate_subsystem
 
+"""
+    FMIComponent(version; fmu, type, name, kwargs...)
+
+Construct a ModelingToolkit component from a Functional Mock-up Unit (FMU).
+
+# Arguments
+
+- `version`: a `Val` selecting the FMI major version, for example `Val(2)`.
+
+# Keyword Arguments
+
+- `fmu`: the loaded FMU object.
+- `type`: the FMI interface type, such as `:ME` for model exchange.
+- `name`: the component name used by `@named`.
+- `kwargs...`: additional FMI extension options.
+
+# Returns
+
+An `AbstractSystem` component representing the FMU interface.
+
+# Examples
+
+```julia
+using ModelingToolkit
+
+@named model = ModelingToolkit.FMIComponent(Val(2); fmu, type = :ME)
+```
+
+!!! note
+    `FMIComponent` is implemented by the FMI extension. Load the FMI dependency before
+    constructing FMU-backed components.
+"""
 function FMIComponent end
 
 @public linearize_symbolic, reorder_unknowns
diff --git a/src/linearization.jl b/src/linearization.jl
index ec130a4dee..b8b53e44c3 100644
--- a/src/linearization.jl
+++ b/src/linearization.jl
@@ -707,14 +707,42 @@ function CommonSolve.solve(prob::LinearizationProblem; allow_input_derivatives =
 end
 
 """
-    (; A, B, C, D), simplified_sys = linearize_symbolic(sys::AbstractSystem, inputs, outputs; simplify = false, allow_input_derivatives = false, kwargs...)
+    matrices, simplified_sys = linearize_symbolic(sys::AbstractSystem, inputs, outputs;
+        simplify = false,
+        allow_input_derivatives = false,
+        eval_expression = false,
+        eval_module = @__MODULE__,
+        split = true,
+        kwargs...)
 
-Similar to [`linearize`](@ref), but returns symbolic matrices `A,B,C,D` rather than numeric. While `linearize` uses ForwardDiff to perform the linearization, this function uses `Symbolics.jacobian`.
+Symbolically linearize `sys` between `inputs` and `outputs`.
 
-See [`linearize`](@ref) for a description of the arguments.
+Unlike [`linearize`](@ref), this uses `Symbolics.jacobian` and returns symbolic
+matrices instead of numerical matrices evaluated at an operating point.
+
+# Arguments
+
+- `sys`: the system to linearize. This function calls [`mtkcompile`](@ref) internally.
+- `inputs`: input variables used as the columns of the `B` and `D` matrices.
+- `outputs`: output variables used as the rows of the `C` and `D` matrices.
+
+# Keyword Arguments
+
+- `simplify`: whether to run additional symbolic simplification during compilation.
+- `allow_input_derivatives`: whether differentiated inputs may appear in the
+  linearized equations. If `true`, `B` and `D` include extra columns for those derivatives.
+- `eval_expression`: whether generated expressions are evaluated directly instead of
+  using runtime-generated functions.
+- `eval_module`: the module used when `eval_expression = true`.
+- `split`: whether generated functions use a tuple of parameters or splatted parameters.
+- `kwargs...`: additional keyword arguments forwarded to [`mtkcompile`](@ref).
+
+# Returns
+
+A pair `(matrices, simplified_sys)`. `matrices` is a `NamedTuple` containing `A`,
+`B`, `C`, and `D`, plus the symbolic Jacobian blocks `f_x`, `f_z`, `g_x`, `g_z`,
+`f_u`, `g_u`, `h_x`, `h_z`, and `h_u` of
 
-# Extended help
-The named tuple returned as the first argument additionally contains the jacobians `f_x, f_z, g_x, g_z, f_u, g_u, h_x, h_z, h_u` of
 ```math
 \\begin{aligned}
 ẋ &= f(x, z, u) \\\\
@@ -722,7 +750,27 @@ ẋ &= f(x, z, u) \\\\
 y &= h(x, z, u)
 \\end{aligned}
 ```
-where `x` are differential unknown variables, `z` algebraic variables, `u` inputs and `y` outputs.
+
+Here `x` are differential unknowns, `z` are algebraic unknowns, `u` are inputs,
+and `y` are outputs.
+
+# Examples
+
+```julia
+using ModelingToolkit
+using ModelingToolkit: t_nounits as t, D_nounits as D
+
+@variables x(t) = 1.0 u(t) = 0.0 y(t) = 0.0
+eqs = [D(x) ~ -x + u,
+       y ~ x]
+@named sys = System(eqs, t)
+
+matrices, simplified_sys = ModelingToolkit.linearize_symbolic(sys, [u], [y])
+matrices.A
+```
+
+See also [`linearize`](@ref), [`linearization_function`](@ref), and
+[`reorder_unknowns`](@ref).
 """
 function linearize_symbolic(
         sys::AbstractSystem, inputs,
@@ -1111,7 +1159,21 @@ end
 """
     (; Ã, B̃, C̃, D̃) = similarity_transform(sys, T; unitary=false)
 
-Perform a similarity transform `T : Tx̃ = x` on linear system represented by matrices in NamedTuple `sys` such that
+Transform the state coordinates of a linear state-space model.
+
+# Arguments
+
+- `sys`: a `NamedTuple` with matrices `A`, `B`, `C`, and `D`, as returned by
+  [`linearize`](@ref) or [`linearize_symbolic`](@ref).
+- `T`: the coordinate transformation matrix where `T * x̃ = x`.
+
+# Keyword Arguments
+
+- `unitary`: if `true`, use the adjoint `T'` instead of factoring and solving with `T`.
+
+# Returns
+
+A `NamedTuple` `(; A, B, C, D)` containing
 
 ```
 Ã = T⁻¹AT
@@ -1120,7 +1182,15 @@ C̃ = CT
 D̃ = D
 ```
 
-If `unitary=true`, `T` is assumed unitary and the matrix adjoint is used instead of the inverse.
+# Examples
+
+```julia
+sys = (; A = [-1.0 0.0; 0.0 -2.0], B = [1.0; 0.0;;],
+       C = [0.0 1.0], D = zeros(1, 1))
+T = [0.0 1.0; 1.0 0.0]
+
+transformed = ModelingToolkit.similarity_transform(sys, T; unitary = true)
+```
 """
 function similarity_transform(sys::NamedTuple, T; unitary = false)
     if unitary
@@ -1139,16 +1209,28 @@ end
 """
     reorder_unknowns(sys::NamedTuple, old, new)
 
-Permute the state representation of `sys` obtained from [`linearize`](@ref) so that the state unknown is changed from `old` to `new`
-Example:
+Permute the state ordering of a linearized system.
 
-```
+# Arguments
+
+- `sys`: a `NamedTuple` with matrices `A`, `B`, `C`, and `D`, as returned by
+  [`linearize`](@ref) or [`linearize_symbolic`](@ref).
+- `old`: the current state order, typically `unknowns(simplified_sys)`.
+- `new`: the desired state order. It must contain the same entries as `old`.
+
+# Returns
+
+A `NamedTuple` `(; A, B, C, D)` whose state-space matrices use the order `new`.
+
+# Examples
+
+```julia
 lsys, ssys = linearize(pid, [reference.u, measurement.u], [ctr_output.u])
 desired_order = [int.x, der.x] # Unknowns that are present in unknowns(ssys)
 lsys = ModelingToolkit.reorder_unknowns(lsys, unknowns(ssys), desired_order)
 ```
 
-See also [`ModelingToolkit.similarity_transform`](@ref)
+See also [`similarity_transform`](@ref).
 """
 function reorder_unknowns(sys::NamedTuple, old, new)
     nx = length(old)
diff --git a/src/problems/docs.jl b/src/problems/docs.jl
index a099793b1c..5e9739bddd 100644
--- a/src/problems/docs.jl
+++ b/src/problems/docs.jl
@@ -26,42 +26,3 @@ Note that all three of `stiff_linear`, `stiff_quadratic`, `stiff_nonlinear` cann
 two of `A`, `B`, `C` returned from [`calculate_semiquadratic_form`](@ref) must be
 non-`nothing`. In other words, both of the functions in the split form must be non-empty.
 """
-
-for (mod, prob, func, istd, kws) in [
-        (SciMLBase, :SCCNonlinearProblem, NonlinearFunction, false, (; init = false)),
-        (
-            ModelingToolkit,
-            :SemilinearODEProblem,
-            :SemilinearODEFunction,
-            true,
-            (;
-                extra_body = SEMILINEAR_EXTRA_BODY, extra_kwargs = SEMILINEAR_A_B_C_KWARGS,
-                extra_kwargs_desc = SEMILINEAR_A_B_C_CONSTRAINT,
-            ),
-        ),
-    ]
-    kwexpr = Expr(:parameters)
-    for (k, v) in pairs(kws)
-        push!(kwexpr.args, Expr(:kw, k, v))
-    end
-    @eval @doc MTKBase.problem_docstring($kwexpr, $mod.$prob, $func, $istd) $mod.$prob
-end
-
-for (mod, func, istd, optionals, kws) in [
-        (
-            ModelingToolkit,
-            :SemilinearODEFunction,
-            true,
-            [:jac],
-            (;
-                extra_body = SEMILINEAR_EXTRA_BODY, extra_kwargs = SEMILINEAR_A_B_C_KWARGS,
-                extra_kwargs_desc = SEMILINEAR_A_B_C_CONSTRAINT,
-            ),
-        ),
-    ]
-    kwexpr = Expr(:parameters)
-    for (k, v) in pairs(kws)
-        push!(kwexpr.args, Expr(:kw, k, v))
-    end
-    @eval @doc MTKBase.function_docstring($kwexpr, $mod.$func, $istd, $optionals) $mod.$func
-end
diff --git a/src/problems/sccnonlinearproblem.jl b/src/problems/sccnonlinearproblem.jl
index 5b555deddc..9769ef6cab 100644
--- a/src/problems/sccnonlinearproblem.jl
+++ b/src/problems/sccnonlinearproblem.jl
@@ -503,6 +503,7 @@ function SCCNonlinearFunction{iip}(
     return NonlinearFunction{iip}(f; sys = subsys)
 end
 
+"""$(MTKBase.problem_docstring(SciMLBase.SCCNonlinearProblem, NonlinearFunction, false; init = false))"""
 function SciMLBase.SCCNonlinearProblem(sys::System, op; kwargs...)
     return SCCNonlinearProblem{true}(sys, op; kwargs...)
 end
diff --git a/src/problems/semilinearodeproblem.jl b/src/problems/semilinearodeproblem.jl
index 10ae4415d7..05db6572dd 100644
--- a/src/problems/semilinearodeproblem.jl
+++ b/src/problems/semilinearodeproblem.jl
@@ -1,3 +1,4 @@
+"""$(MTKBase.function_docstring(SemilinearODEFunction, true, [:jac]; extra_body = SEMILINEAR_EXTRA_BODY, extra_kwargs = SEMILINEAR_A_B_C_KWARGS, extra_kwargs_desc = SEMILINEAR_A_B_C_CONSTRAINT))"""
 @fallback_iip_specialize function SemilinearODEFunction{iip, specialize}(
         sys::System; u0 = nothing, p = nothing, t = nothing,
         semiquadratic_form = nothing,
@@ -88,6 +89,7 @@ function SemilinearODEFunction{iip, specialize}(
     )
 end
 
+"""$(MTKBase.problem_docstring(SemilinearODEProblem, SemilinearODEFunction, true; extra_body = SEMILINEAR_EXTRA_BODY, extra_kwargs = SEMILINEAR_A_B_C_KWARGS, extra_kwargs_desc = SEMILINEAR_A_B_C_CONSTRAINT))"""
 @fallback_iip_specialize function SemilinearODEProblem{iip, spec}(
         sys::System, op, tspan; check_compatibility = true, u0_eltype = nothing,
         expression = Val{false}, callback = nothing, sparse = false,
diff --git a/src/structural_transformation/StructuralTransformations.jl b/src/structural_transformation/StructuralTransformations.jl
index 4af214d4a5..2ecace3fd5 100644
--- a/src/structural_transformation/StructuralTransformations.jl
+++ b/src/structural_transformation/StructuralTransformations.jl
@@ -68,9 +68,21 @@ include("utils.jl")
 include("pantelides.jl")
 
 """
-    $TYPEDSIGNATURES
+    tearing_substitution(sys::AbstractSystem; kwargs...)
 
-Apply tearing substitutions to `sys` and clear its schedule.
+Replace the equations of `sys` with its fully substituted equations.
+
+This is a structural-transformation helper used by simplification passes. End-user code
+should usually call [`mtkcompile`](@ref).
+
+# Arguments
+
+- `sys`: system whose equations should be substituted.
+- `kwargs...`: keyword arguments forwarded to [`full_equations`](@ref).
+
+# Returns
+
+A copy of `sys` with substituted equations and no cached schedule.
 """
 function tearing_substitution(sys::AbstractSystem; kwargs...)
     neweqs = full_equations(sys::AbstractSystem; kwargs...)
diff --git a/src/structural_transformation/pantelides.jl b/src/structural_transformation/pantelides.jl
index c0f0be3f6b..aa477b2667 100644
--- a/src/structural_transformation/pantelides.jl
+++ b/src/structural_transformation/pantelides.jl
@@ -5,9 +5,19 @@
 const NOTHING_EQ = nothing ~ nothing
 
 """
-    $TYPEDSIGNATURES
+    pantelides_reassemble(state::TearingState, var_eq_matching)
 
-Reassemble a system from Pantelides structural state and a variable-equation matching.
+Reassemble a `System` after Pantelides index reduction.
+
+# Arguments
+
+- `state`: tearing state containing the original system and derivative graphs.
+- `var_eq_matching`: variable-equation matching returned by the Pantelides pass.
+
+# Returns
+
+A system whose equations and unknowns include the differentiated equations selected by
+the Pantelides algorithm.
 """
 function pantelides_reassemble(state::TearingState, var_eq_matching)
     fullvars = state.fullvars
diff --git a/src/structural_transformation/utils.jl b/src/structural_transformation/utils.jl
index 87c0019ef7..7037db16aa 100644
--- a/src/structural_transformation/utils.jl
+++ b/src/structural_transformation/utils.jl
@@ -27,9 +27,17 @@ end
 ### Structural and symbolic utilities
 ###
 """
-    $TYPEDSIGNATURES
+    highest_order_variable_mask(ts)
 
-Return a predicate selecting variables with no differentiated successor in the tearing state.
+Return a predicate selecting highest-order variables in a tearing state.
+
+# Arguments
+
+- `ts`: tearing state whose derivative graph should be inspected.
+
+# Returns
+
+A predicate `f(v)::Bool` over variable indices.
 """
 function highest_order_variable_mask(ts)
     return let v2d = ts.structure.var_to_diff
@@ -38,9 +46,17 @@ function highest_order_variable_mask(ts)
 end
 
 """
-    $TYPEDSIGNATURES
+    lowest_order_variable_mask(ts)
+
+Return a predicate selecting lowest-order variables in a tearing state.
+
+# Arguments
 
-Return a predicate selecting variables with no lower-order predecessor in the tearing state.
+- `ts`: tearing state whose derivative graph should be inspected.
+
+# Returns
+
+A predicate `f(v)::Bool` over variable indices.
 """
 function lowest_order_variable_mask(ts)
     return let v2d = ts.structure.var_to_diff
@@ -49,9 +65,19 @@ function lowest_order_variable_mask(ts)
 end
 
 """
-    $TYPEDSIGNATURES
+    but_ordered_incidence(ts::TearingState, varmask = highest_order_variable_mask(ts))
+
+Construct the block upper triangular ordered incidence matrix for `ts`.
+
+# Arguments
+
+- `ts`: tearing state to analyze.
+- `varmask`: predicate selecting which variable indices participate in the ordering.
+
+# Returns
 
-Return a block upper triangular incidence matrix and block boundaries for `ts`.
+A pair `(matrix, block_boundaries)` containing the ordered incidence matrix and the
+starting row/column boundaries for each block.
 """
 function but_ordered_incidence(ts::TearingState, varmask = highest_order_variable_mask(ts))
     graph = complete(ts.structure.graph)
diff --git a/src/systems/alias_elimination.jl b/src/systems/alias_elimination.jl
index 00a7ca301b..9446972fdb 100644
--- a/src/systems/alias_elimination.jl
+++ b/src/systems/alias_elimination.jl
@@ -2,9 +2,20 @@ using SymbolicUtils: Rewriters
 using Graphs.Experimental.Traversals
 
 """
-    $TYPEDSIGNATURES
+    alias_elimination(sys)
+
+Run the alias-elimination pass on `sys`.
+
+This is a lower-level structural simplification pass used by [`mtkcompile`](@ref).
+End-user code should usually call `mtkcompile` instead.
+
+# Arguments
+
+- `sys`: system to transform.
+
+# Returns
 
-Return a system with perfect aliases eliminated.
+The alias-eliminated system.
 """
 alias_elimination(sys) = alias_elimination!(TearingState(sys))[1]
 
diff --git a/src/systems/analysis_points.jl b/src/systems/analysis_points.jl
index dcfa29a8e5..f949bae8cb 100644
--- a/src/systems/analysis_points.jl
+++ b/src/systems/analysis_points.jl
@@ -116,20 +116,93 @@ function get_looptransfer_function(sys::AbstractSystem, aps; kwargs...)
     return get_linear_analysis_function(sys, LoopTransferTransform, aps; kwargs...)
 end
 
-for f in [:get_sensitivity, :get_comp_sensitivity, :get_looptransfer]
-    utility_fun = Symbol(f, :_function)
-    @eval function $f(
-            sys, ap, args...; loop_openings = [], system_modifier = identity,
-            allow_input_derivatives = true, op = Dict{SymbolicT, SymbolicT}(), t = 0.0,
-            kwargs...
-        )
-        lin_fun,
-            ssys = $(utility_fun)(
-            sys, ap, args...; loop_openings, system_modifier, op, kwargs...
-        )
-        mats, extras = ModelingToolkit.linearize(ssys, lin_fun; op, allow_input_derivatives, t)
-        return mats, ssys, extras
-    end
+"""
+    get_sensitivity(sys, ap; loop_openings = [], system_modifier = identity, op = Dict(), kwargs...)
+
+Compute the sensitivity transfer matrices at analysis point `ap` and return
+`(matrices, simplified_system, extras)`. This introduces an infinitesimal perturbation
+at the input of `ap`, linearizes the system, and computes the transfer function from the
+perturbation to the output of `ap`.
+
+# Arguments
+
+- `sys`: system containing the analysis point.
+- `ap`: an [`AnalysisPoint`](@ref) or its symbolic name.
+
+# Keyword Arguments
+
+- `loop_openings`: connections to open before linearization.
+- `system_modifier`: transformation applied before linearization.
+- `op`: operating-point values passed to [`linearize`](@ref).
+- `allow_input_derivatives`: allow derivatives of input variables in the linearization.
+- `t`: time at which to evaluate the linearization.
+- `kwargs...`: forwarded to [`get_sensitivity_function`](@ref).
+
+See also [`get_comp_sensitivity`](@ref) and [`get_looptransfer`](@ref).
+"""
+function get_sensitivity(
+        sys, ap, args...; loop_openings = [], system_modifier = identity,
+        allow_input_derivatives = true, op = Dict{SymbolicT, SymbolicT}(), t = 0.0,
+        kwargs...
+    )
+    lin_fun, ssys = get_sensitivity_function(
+        sys, ap, args...; loop_openings, system_modifier, op, kwargs...
+    )
+    mats, extras = ModelingToolkit.linearize(ssys, lin_fun; op, allow_input_derivatives, t)
+    return mats, ssys, extras
+end
+
+"""
+    get_comp_sensitivity(sys, ap; loop_openings = [], system_modifier = identity, op = Dict(), kwargs...)
+
+Compute the complementary-sensitivity transfer matrices at analysis point `ap` and return
+`(matrices, simplified_system, extras)`. This introduces an infinitesimal perturbation
+at the output of `ap`, linearizes the system, and computes the transfer function from the
+perturbation to the input of `ap`.
+
+Keyword arguments match [`get_sensitivity`](@ref), with `kwargs...` forwarded to
+[`get_comp_sensitivity_function`](@ref).
+
+See also [`get_sensitivity`](@ref) and [`get_looptransfer`](@ref).
+"""
+function get_comp_sensitivity(
+        sys, ap, args...; loop_openings = [], system_modifier = identity,
+        allow_input_derivatives = true, op = Dict{SymbolicT, SymbolicT}(), t = 0.0,
+        kwargs...
+    )
+    lin_fun, ssys = get_comp_sensitivity_function(
+        sys, ap, args...; loop_openings, system_modifier, op, kwargs...
+    )
+    mats, extras = ModelingToolkit.linearize(ssys, lin_fun; op, allow_input_derivatives, t)
+    return mats, ssys, extras
+end
+
+"""
+    get_looptransfer(sys, ap; loop_openings = [], system_modifier = identity, op = Dict(), kwargs...)
+
+Compute the loop-transfer matrices at analysis point `ap` and return
+`(matrices, simplified_system, extras)`. The transfer is from `ap.out` to `ap.in`.
+
+Keyword arguments match [`get_sensitivity`](@ref), with `kwargs...` forwarded to
+[`get_looptransfer_function`](@ref).
+
+!!! info "Negative feedback"
+
+    The computed loop transfer includes negative feedback. Negate the result when using an
+    analysis tool that expects a loop-transfer function without the negative gain.
+
+See also [`get_sensitivity`](@ref), [`get_comp_sensitivity`](@ref), and [`open_loop`](@ref).
+"""
+function get_looptransfer(
+        sys, ap, args...; loop_openings = [], system_modifier = identity,
+        allow_input_derivatives = true, op = Dict{SymbolicT, SymbolicT}(), t = 0.0,
+        kwargs...
+    )
+    lin_fun, ssys = get_looptransfer_function(
+        sys, ap, args...; loop_openings, system_modifier, op, kwargs...
+    )
+    mats, extras = ModelingToolkit.linearize(ssys, lin_fun; op, allow_input_derivatives, t)
+    return mats, ssys, extras
 end
 
 """
@@ -476,46 +549,3 @@ function isolate_subsystem(
 
     return _reconstruct!(sys, Symbol[], all_clock_subs), input_vars, output_vars
 end
-
-@doc """
-    get_sensitivity(sys, ap::AnalysisPoint; kwargs)
-    get_sensitivity(sys, ap_name::Symbol; kwargs)
-
-Compute the sensitivity function in analysis point `ap`. The sensitivity function is obtained by introducing an infinitesimal perturbation `d` at the input of `ap`, linearizing the system and computing the transfer function between `d` and the output of `ap`.
-
-# Arguments:
-
-  - `kwargs`: Are sent to `ModelingToolkit.linearize`
-
-See also [`get_comp_sensitivity`](@ref), [`get_looptransfer`](@ref).
-""" get_sensitivity
-
-@doc """
-    get_comp_sensitivity(sys, ap::AnalysisPoint; kwargs)
-    get_comp_sensitivity(sys, ap_name::Symbol; kwargs)
-
-Compute the complementary sensitivity function in analysis point `ap`. The complementary sensitivity function is obtained by introducing an infinitesimal perturbation `d` at the output of `ap`, linearizing the system and computing the transfer function between `d` and the input of `ap`.
-
-# Arguments:
-
-  - `kwargs`: Are sent to `ModelingToolkit.linearize`
-
-See also [`get_sensitivity`](@ref), [`get_looptransfer`](@ref).
-""" get_comp_sensitivity
-
-@doc """
-    get_looptransfer(sys, ap::AnalysisPoint; kwargs)
-    get_looptransfer(sys, ap_name::Symbol; kwargs)
-
-Compute the (linearized) loop-transfer function in analysis point `ap`, from `ap.out` to `ap.in`.
-
-!!! info "Negative feedback"
-
-    Feedback loops often use negative feedback, and the computed loop-transfer function will in this case have the negative feedback included. Standard analysis tools often assume a loop-transfer function without the negative gain built in, and the result of this function may thus need negation before use.
-
-# Arguments:
-
-  - `kwargs`: Are sent to `ModelingToolkit.linearize`
-
-See also [`get_sensitivity`](@ref), [`get_comp_sensitivity`](@ref), [`open_loop`](@ref).
-""" get_looptransfer
diff --git a/src/systems/substitute_component.jl b/src/systems/substitute_component.jl
index 122cd75a8a..0841875044 100644
--- a/src/systems/substitute_component.jl
+++ b/src/systems/substitute_component.jl
@@ -151,6 +151,14 @@ Given a hierarchical system `sys` and a rule `lhs => rhs`, replace the subsystem
 obtained via `sys.inner.component`). The `rhs` must be valid as per the following
 conditions:
 
+# Arguments
+
+- `sys`: the system whose nested subsystem should be replaced.
+- `rule`: a pair `lhs => rhs`, where `lhs` is the namespaced subsystem to replace and
+  `rhs` is the unnamespaced replacement system.
+
+# Replacement Rules
+
 1. `rhs` must not be namespaced.
 2. The name of `rhs` must be the same as the unnamespaced name of `lhs`.
 3. Neither one of `lhs` or `rhs` can be marked as complete.
@@ -165,6 +173,29 @@ conditions:
    be a connector of the same type.
 
 `sys` also cannot be marked as complete.
+
+# Returns
+
+A new system with the selected subsystem replaced and connection equations rebuilt.
+
+# Examples
+
+```julia
+using ModelingToolkit
+using ModelingToolkit: t_nounits as t, D_nounits as D
+
+function Block(k; name)
+    @variables x(t) = 0
+    @parameters p = k
+    System([D(x) ~ -p * x], t; name)
+end
+
+@named old_block = Block(1)
+@named new_block = Block(2)
+@named sys = System(Equation[], t; systems = [old_block])
+
+updated = substitute_component(sys, sys.old_block => new_block)
+```
 """
 function substitute_component(sys::T, rule::Pair{T, T}) where {T <: AbstractSystem}
     iscomplete(sys) &&
diff --git a/src/systems/systems.jl b/src/systems/systems.jl
index f9ffea4f78..c31d6cb429 100644
--- a/src/systems/systems.jl
+++ b/src/systems/systems.jl
@@ -1,25 +1,3 @@
-@doc """
-    function mtkcompile(sys::System; kwargs...)
-
-Compile the given system into a form that ModelingToolkit can generate code for. Also
-performs a variety of symbolic-numeric enhancements. For ODEs, this includes processes
-such as order reduction, index reduction, alias elimination and tearing. A subset of the
-unknowns of the system may be eliminated as observables, eliminating the need for the
-numerical solver to solve for these variables.
-
-Does not rely on metadata to identify variables/parameters/brownians. Instead, queries
-the system for which symbolic quantities belong to which category. Any variables not
-present in the equations of the system will be removed in this process.
-
-# Keyword Arguments
-
-+ When `simplify=true`, the `simplify` function will be applied during the tearing process.
-+ `allow_symbolic=false`, `allow_parameter=true`, and `conservative=false` limit the coefficient types during tearing. In particular, `conservative=true` limits tearing to only solve for trivial linear systems where the coefficient has the absolute value of ``1``.
-+ `fully_determined=true` controls whether or not an error will be thrown if the number of equations don't match the number of inputs, outputs, and equations.
-+ `inputs`, `outputs` and `disturbance_inputs` are passed as keyword arguments.` All inputs` get converted to parameters and are allowed to be unconnected, allowing models where `n_unknowns = n_equations - n_inputs`.
-+ `sort_eqs=true` controls whether equations are sorted lexicographically before simplification or not.
-""" mtkcompile
-
 function MTKBase.__mtkcompile(
         sys::System;
         inputs::OrderedSet{SymbolicT} = OrderedSet{SymbolicT}(),
@@ -161,11 +139,37 @@ Given a system that has been simplified via `mtkcompile`, return a `Dict` mappin
 variables of the system to equations that are used to solve for them. This includes
 observed variables.
 
+# Arguments
+
+- `sys`: a system returned by [`mtkcompile`](@ref) or another simplification path that
+  records a tearing state.
+
 # Keyword Arguments
 
 - `rename_dummy_derivatives`: Whether to rename dummy derivative variable keys into their
   `Differential` forms. For example, this would turn the key `yˍt(t)` into
   `Differential(t)(y(t))`.
+
+# Returns
+
+A `Dict` mapping unknown or observed symbolic variables to the equations used to solve
+for them.
+
+# Examples
+
+```julia
+using ModelingToolkit
+using ModelingToolkit: t_nounits as t, D_nounits as D
+
+@variables x(t) = 1 y(t) = 0
+eqs = [D(x) ~ -x,
+       y ~ x + 1]
+@named sys = System(eqs, t)
+simplified = mtkcompile(sys)
+
+mapping = map_variables_to_equations(simplified)
+mapping[y]
+```
 """
 function map_variables_to_equations(sys::AbstractSystem; rename_dummy_derivatives = true)
     if !has_tearing_state(sys)
diff --git a/src/systems/systemstructure.jl b/src/systems/systemstructure.jl
index 60b9e0d8fb..c41a5373f7 100644
--- a/src/systems/systemstructure.jl
+++ b/src/systems/systemstructure.jl
@@ -98,10 +98,31 @@ function inputs_to_parameters!(state::TearingState, inputsyms::OrderedSet{Symbol
 end
 
 """
-    $TYPEDSIGNATURES
+    mtkcompile!(state::TearingState; kwargs...)
+
+Mutating structural simplification entry point for an existing tearing state.
+
+This is developer-facing API used by ModelingToolkit internals and extension packages that
+already have a `TearingState`. User code should normally call [`mtkcompile`](@ref) on a
+`System` instead.
+
+# Arguments
+
+- `state`: tearing state to simplify in place.
+
+# Keyword Arguments
+
+- `check_consistency`: whether to check the transformed system for structural consistency.
+- `fully_determined`: whether the transformed system is expected to have a square
+  equation/unknown structure.
+- `inputs`: variables to treat as external inputs.
+- `outputs`: variables to treat as requested outputs.
+- `disturbance_inputs`: input variables that should be treated as disturbances.
+- `kwargs...`: additional simplification options forwarded to the internal compiler.
+
+# Returns
 
-Compile the system stored in `state` in place, updating its tearing state and returning the
-compiled [`System`](@ref).
+The simplified `System`.
 """
 function mtkcompile!(
         state::TearingState;
diff --git a/test/downstream/Project.toml b/test/downstream/Project.toml
index 0bfba05dfa..a22cf7e1e1 100644
--- a/test/downstream/Project.toml
+++ b/test/downstream/Project.toml
@@ -19,4 +19,4 @@ SymbolicIndexingInterface = "2efcf032-c050-4f8e-a9bb-153293bab1f5"
 [compat]
 ModelingToolkitStandardLibrary = "2.19"
 SafeTestsets = "0.1, 1"
-SciMLTesting = "1, 2.1"
+SciMLTesting = "2.8"
diff --git a/test/fmi/Project.toml b/test/fmi/Project.toml
index ae5bc6dce9..d5840ab563 100644
--- a/test/fmi/Project.toml
+++ b/test/fmi/Project.toml
@@ -9,4 +9,4 @@ SciMLTesting = "09d9d899-5365-40a9-917a-5f67fddea283"
 [compat]
 FMI = "0.14"
 SafeTestsets = "0.1, 1"
-SciMLTesting = "1, 2.1"
+SciMLTesting = "2.8"
diff --git a/test/qa/Project.toml b/test/qa/Project.toml
index e8786faa58..8d133366aa 100644
--- a/test/qa/Project.toml
+++ b/test/qa/Project.toml
@@ -27,7 +27,7 @@ OrdinaryDiffEqBDF = "1, 2"
 OrdinaryDiffEqDefault = "1.2, 2"
 OrdinaryDiffEqRosenbrock = "1, 2"
 SafeTestsets = "0.1, 1"
-SciMLTesting = "1, 2.1"
+SciMLTesting = "2.8"
 StateSelection = "1.10.5"
 Test = "1"
 julia = "1.10"