diff --git a/docs/Project.toml b/docs/Project.toml
index 314a1d7f84..746d6448b6 100644
--- a/docs/Project.toml
+++ b/docs/Project.toml
@@ -3,6 +3,7 @@ Attractors = "f3fd9213-ca85-4dba-9dfd-7fc91308fec7"
 BenchmarkTools = "6e4b80f9-dd63-53aa-95a3-0cdb28fa8baf"
 BifurcationKit = "0f109fa4-8a5d-4b75-95aa-f515264e7665"
 CairoMakie = "13f3f980-e62b-5c42-98c6-ff1f3baf88f0"
+CommonSolve = "38540f10-b2f7-11e9-35d8-d573e4eb0ff2"
 ControlSystemsBase = "aaaaaaaa-a6ca-5380-bf3e-84a91bcd477e"
 DataInterpolations = "82cc6244-b520-54b8-b5a6-8a565e85f1d0"
 DiffEqDevTools = "f3b72e0c-5b89-59e1-b016-84e28bfd966d"
@@ -13,6 +14,7 @@ FMI = "14a09403-18e3-468f-ad8a-74f8dda2d9ac"
 FMIZoo = "724179cf-c260-40a9-bd27-cccc6fe2f195"
 InfiniteOpt = "20393b10-9daf-11e9-18c9-8db751c92c57"
 Ipopt = "b6b21f68-93f8-5de0-b562-5493be1d77c9"
+JumpProcesses = "ccbc3e58-028d-4f4c-8cd5-9ae44345cda5"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 ModelingToolkit = "961ee093-0014-501f-94e3-6117800e7a78"
 ModelingToolkitStandardLibrary = "16a59e39-deab-5bd0-87e4-056b12336739"
diff --git a/docs/make.jl b/docs/make.jl
index 36a27bf598..248b9e35bc 100644
--- a/docs/make.jl
+++ b/docs/make.jl
@@ -1,5 +1,9 @@
 using Documenter, ModelingToolkit
 using ModelingToolkit: SciMLBase
+# To load docstring from extension
+import FMI, CommonSolve, JumpProcesses
+
+MTKFMIExt = Base.get_extension(ModelingToolkit, :MTKFMIExt)
 
 # Make sure that plots don't throw a bunch of warnings / errors!
 ENV["GKSwstype"] = "100"
@@ -22,7 +26,7 @@ mathengine = MathJax3(Dict(:loader => Dict("load" => ["[tex]/require", "[tex]/ma
 
 makedocs(sitename = "ModelingToolkit.jl",
     authors = "Chris Rackauckas",
-    modules = [ModelingToolkit],
+    modules = [ModelingToolkit, MTKFMIExt],
     clean = true, doctest = false, linkcheck = true,
     warnonly = [:docs_block, :missing_docs, :cross_references],
     linkcheck_ignore = [
@@ -35,7 +39,9 @@ makedocs(sitename = "ModelingToolkit.jl",
         assets = ["assets/favicon.ico"],
         mathengine,
         canonical = "https://docs.sciml.ai/ModelingToolkit/stable/",
-        prettyurls = (get(ENV, "CI", nothing) == "true")),
+        prettyurls = (get(ENV, "CI", nothing) == "true"),
+        # This page gets especially big with all the problem docstrings
+        size_threshold_ignore = ["API/problems.md"]),
     pages = pages)
 
 deploydocs(repo = "github.com/SciML/ModelingToolkit.jl.git";
diff --git a/docs/src/API/System.md b/docs/src/API/System.md
index 33a4d064ab..379e7fa19b 100644
--- a/docs/src/API/System.md
+++ b/docs/src/API/System.md
@@ -32,7 +32,7 @@ or analysis points of the hierarchical system.
 ModelingToolkit.has_eqs
 ModelingToolkit.get_eqs
 equations
-equations_toplevel
+ModelingToolkit.equations_toplevel
 full_equations
 ModelingToolkit.has_noise_eqs
 ModelingToolkit.get_noise_eqs
@@ -44,17 +44,17 @@ ModelingToolkit.get_constraints
 constraints
 ModelingToolkit.has_costs
 ModelingToolkit.get_costs
-costs
+cost
 ModelingToolkit.has_consolidate
 ModelingToolkit.get_consolidate
 ModelingToolkit.has_unknowns
 ModelingToolkit.get_unknowns
 unknowns
-unknowns_toplevel
+ModelingToolkit.unknowns_toplevel
 ModelingToolkit.has_ps
 ModelingToolkit.get_ps
 parameters
-parameters_toplevel
+ModelingToolkit.parameters_toplevel
 tunable_parameters
 ModelingToolkit.has_brownians
 ModelingToolkit.get_brownians
@@ -77,19 +77,20 @@ defaults
 ModelingToolkit.has_guesses
 ModelingToolkit.get_guesses
 guesses
+ModelingToolkit.get_systems
 ModelingToolkit.has_initialization_eqs
 ModelingToolkit.get_initialization_eqs
 initialization_equations
 ModelingToolkit.has_continuous_events
 ModelingToolkit.get_continuous_events
 continuous_events
-continuous_events_toplevel
+ModelingToolkit.continuous_events_toplevel
 ModelingToolkit.has_discrete_events
 ModelingToolkit.get_discrete_events
-discrete_events_toplevel
+ModelingToolkit.discrete_events_toplevel
 ModelingToolkit.has_assertions
 ModelingToolkit.get_assertions
-assertions
+ModelingToolkit.assertions
 ModelingToolkit.has_metadata
 ModelingToolkit.get_metadata
 SymbolicUtils.getmetadata(::ModelingToolkit.AbstractSystem, ::DataType, ::Any)
@@ -144,8 +145,8 @@ has_diff_equations
 has_alg_equations
 diff_equations
 alg_equations
-is_alg_equation
-is_diff_equation
+ModelingToolkit.is_alg_equation
+ModelingToolkit.is_diff_equation
 ```
 
 ## String parsing
@@ -167,6 +168,18 @@ ModelingToolkit.dump_parameters
 ModelingToolkit.dump_variable_metadata
 ```
 
+## Inputs and outputs
+
+```@docs
+ModelingToolkit.inputs
+ModelingToolkit.outputs
+ModelingToolkit.bound_inputs
+ModelingToolkit.unbound_inputs
+ModelingToolkit.bound_outputs
+ModelingToolkit.unbound_outputs
+ModelingToolkit.is_bound
+```
+
 ## Debugging utilities
 
 ```@docs
diff --git a/docs/src/API/codegen.md b/docs/src/API/codegen.md
index b3a9c01e58..cd76d34522 100644
--- a/docs/src/API/codegen.md
+++ b/docs/src/API/codegen.md
@@ -8,7 +8,6 @@ ModelingToolkit.generate_rhs
 ModelingToolkit.generate_diffusion_function
 ModelingToolkit.generate_jacobian
 ModelingToolkit.generate_tgrad
-ModelingToolkit.generate_hessian
 ModelingToolkit.generate_W
 ModelingToolkit.generate_dae_jacobian
 ModelingToolkit.generate_history
@@ -21,6 +20,7 @@ ModelingToolkit.generate_constraint_jacobian
 ModelingToolkit.generate_constraint_hessian
 ModelingToolkit.generate_control_jacobian
 ModelingToolkit.build_explicit_observed_function
+ModelingToolkit.generate_control_function
 ```
 
 For functions such as jacobian calculation which require symbolic computation, there
@@ -43,3 +43,9 @@ ModelingToolkit.calculate_constraint_jacobian
 ModelingToolkit.calculate_constraint_hessian
 ModelingToolkit.calculate_control_jacobian
 ```
+
+All code generation eventually calls `build_function_wrapper`.
+
+```@docs
+build_function_wrapper
+```
diff --git a/docs/src/API/dynamic_opt.md b/docs/src/API/dynamic_opt.md
index 1a7081154c..b086fca522 100644
--- a/docs/src/API/dynamic_opt.md
+++ b/docs/src/API/dynamic_opt.md
@@ -28,7 +28,7 @@ JuMPCollocation
 InfiniteOptCollocation
 CasADiCollocation
 PyomoCollocation
-solve(::AbstractDynamicOptProblem)
+CommonSolve.solve(::AbstractDynamicOptProblem)
 ```
 
 ### Problem constructors
diff --git a/docs/src/API/model_building.md b/docs/src/API/model_building.md
index 40072bc2c4..a72be8a4c6 100644
--- a/docs/src/API/model_building.md
+++ b/docs/src/API/model_building.md
@@ -81,7 +81,7 @@ Similar to the `stream` and `flow` keyword arguments in the specification, Model
 allows specifying how variables in a connector behave in a connection.
 
 ```@docs
-Equality
+ModelingToolkit.Equality
 Flow
 Stream
 ```
@@ -168,7 +168,6 @@ Symbolic affects are handled using equations as described in the [Events](@ref e
 section of the documentation. User-defined functions can be used via `ImperativeAffect`.
 
 ```@docs
-ModelingToolkit.AffectSystem
 ModelingToolkit.ImperativeAffect
 ```
 
diff --git a/docs/src/API/problems.md b/docs/src/API/problems.md
index 5549224a09..1e3591ea26 100644
--- a/docs/src/API/problems.md
+++ b/docs/src/API/problems.md
@@ -9,46 +9,46 @@ code for a variety of such numerical problems.
 ## Dynamical systems
 
 ```@docs
-ODEFunction(::System, args...)
-ODEProblem(::System, args...)
-DAEFunction(::System, args...)
-DAEProblem(::System, args...)
-SDEFunction(::System, args...)
-SDEProblem(::System, args...)
-DDEFunction(::System, args...)
-DDEProblem(::System, args...)
-SDDEFunction(::System, args...)
-SDDEProblem(::System, args...)
-JumpProblem(::System, args...)
-BVProblem(::System, args...)
-DiscreteProblem(::System, args...)
-ImplicitDiscreteProblem(::System, args...)
+SciMLBase.ODEFunction
+SciMLBase.ODEProblem
+SciMLBase.DAEFunction
+SciMLBase.DAEProblem
+SciMLBase.SDEFunction
+SciMLBase.SDEProblem
+SciMLBase.DDEFunction
+SciMLBase.DDEProblem
+SciMLBase.SDDEFunction
+SciMLBase.SDDEProblem
+JumpProcesses.JumpProblem
+SciMLBase.BVProblem
+SciMLBase.DiscreteProblem
+SciMLBase.ImplicitDiscreteProblem
 ```
 
 ## Nonlinear systems
 
 ```@docs
-NonlinearFunction(::System, args...)
-NonlinearProblem(::System, args...)
-SCCNonlinearProblem(::System, args...)
-NonlinearLeastSquaresProblem(::System, args...)
-SteadyStateProblem(::System, args...)
-IntervalNonlinearFunction(::System, args...)
-IntervalNonlinearProblem(::System, args...)
+SciMLBase.NonlinearFunction
+SciMLBase.NonlinearProblem
+SciMLBase.SCCNonlinearProblem
+SciMLBase.NonlinearLeastSquaresProblem
+SciMLBase.SteadyStateProblem
+SciMLBase.IntervalNonlinearFunction
+SciMLBase.IntervalNonlinearProblem
 ModelingToolkit.HomotopyContinuationProblem
-HomotopyNonlinearFunction(::System, args...)
+SciMLBase.HomotopyNonlinearFunction
 ```
 
 ## Optimization and optimal control
 
 ```@docs
-OptimizationFunction(::System, args...)
-OptimizationProblem(::System, args...)
-ODEInputFunction(::System, args...)
-JuMPDynamicOptProblem(::System, args...)
-InfiniteOptDynamicOptProblem,(::System, args...)
-CasADiDynamicOptProblem(::System, args...)
-DynamicOptSolution
+SciMLBase.OptimizationFunction
+SciMLBase.OptimizationProblem
+SciMLBase.ODEInputFunction
+ModelingToolkit.JuMPDynamicOptProblem
+ModelingToolkit.InfiniteOptDynamicOptProblem
+ModelingToolkit.CasADiDynamicOptProblem
+ModelingToolkit.DynamicOptSolution
 ```
 
 ## The state vector and parameter object
@@ -102,7 +102,7 @@ There are also utilities for manipulating the results of these analyses in a sym
 
 ```@docs
 ModelingToolkit.similarity_transform
-ModelingToolkit.reorder_unknnowns
+ModelingToolkit.reorder_unknowns
 ```
 
 ### Analysis point transformations
diff --git a/docs/src/API/variables.md b/docs/src/API/variables.md
index 04d85e06b9..6f48c9a5f8 100644
--- a/docs/src/API/variables.md
+++ b/docs/src/API/variables.md
@@ -1,4 +1,4 @@
-# Symbolic variables and variable metadata
+# [Symbolic variables and variable metadata](@id symbolic_metadata)
 
 ModelingToolkit uses [Symbolics.jl](https://docs.sciml.ai/Symbolics/stable/) for the symbolic
 manipulation infrastructure. In fact, the `@variables` macro is defined in Symbolics.jl. In
@@ -85,7 +85,7 @@ hasconnect
 getconnect
 ```
 
-```@docs; canonical = false
+```@docs; canonical=false
 Flow
 Stream
 ```
@@ -165,7 +165,7 @@ getguess
 When a system is constructed, the guesses of the involved variables are stored in a `Dict`
 in the system. After this point, the guess metadata of the variable is irrelevant.
 
-```@docs; canonical = false
+```@docs; canonical=false
 guesses
 ```
 
@@ -315,7 +315,7 @@ b = getbounds(sys) # Operating on the system, we get a dict
 
 See also:
 
-```@docs; canonical = false
+```@docs; canonical=false
 tunable_parameters
 ModelingToolkit.dump_unknowns
 ModelingToolkit.dump_parameters
diff --git a/docs/src/basics/AbstractSystem.md b/docs/src/basics/AbstractSystem.md
deleted file mode 100644
index 61b6ef4fff..0000000000
--- a/docs/src/basics/AbstractSystem.md
+++ /dev/null
@@ -1,162 +0,0 @@
-# The AbstractSystem Interface
-
-## Overview
-
-The `AbstractSystem` interface is the core of the system level of ModelingToolkit.jl.
-It establishes a common set of functionality that is used between systems
-representing ODEs, PDEs, SDEs and more, allowing users to have a common framework for
-model manipulation and compilation.
-
-### Subtypes
-
-There are three immediate subtypes of `AbstractSystem`, classified by how many independent variables each type has:
-
-  - `AbstractTimeIndependentSystem`: has no independent variable (e.g.: `NonlinearSystem`)
-  - `AbstractTimeDependentSystem`: has a single independent variable (e.g.: `System`)
-  - `AbstractMultivariateSystem`: may have multiple independent variables (e.g.: `PDESystem`)
-
-## Constructors and Naming
-
-The `AbstractSystem` interface has a consistent method for constructing systems.
-Generally, it follows the order of:
-
- 1. Equations
- 2. Independent Variables
- 3. Dependent Variables (or Unknowns)
- 4. Parameters
-
-All other pieces are handled via keyword arguments. `AbstractSystem`s share the
-same keyword arguments, which are:
-
-  - `system`: This is used for specifying subsystems for hierarchical modeling with
-    reusable components. For more information, see the [components page](@ref components).
-  - Defaults: Keyword arguments like `defaults` are used for specifying default
-    values which are used. If a value is not given at the `SciMLProblem` construction
-    time, its numerical value will be the default.
-
-## Composition and Accessor Functions
-
-Each `AbstractSystem` has lists of variables in context, such as distinguishing
-parameters vs unknowns. In addition, an `AbstractSystem` can also hold other
-`AbstractSystem` types. Direct accessing of the values, such as `sys.unknowns`,
-gives the immediate list, while the accessor functions `unknowns(sys)` gives the
-total set, which includes that of all systems held inside.
-
-The values which are common to all `AbstractSystem`s are:
-
-  - `equations(sys)`: All equations that define the system and its subsystems.
-  - `unknowns(sys)`: All the unknowns in the system and its subsystems.
-  - `parameters(sys)`: All parameters of the system and its subsystems.
-  - `nameof(sys)`: The name of the current-level system.
-  - `get_eqs(sys)`: Equations that define the current-level system.
-  - `get_unknowns(sys)`: Unknowns that are in the current-level system.
-  - `get_ps(sys)`: Parameters that are in the current-level system.
-  - `get_systems(sys)`: Subsystems of the current-level system.
-
-Optionally, a system could have:
-
-  - `observed(sys)`: All observed equations of the system and its subsystems.
-  - `independent_variables(sys)`: The independent variables of a system.
-  - `defaults(sys)`: A `Dict` that maps variables/parameters into their default values for the system and its subsystems.
-  - `get_observed(sys)`: Observed equations of the current-level system.
-  - `get_continuous_events(sys)`: `SymbolicContinuousCallback`s of the current-level system.
-  - `get_defaults(sys)`: A `Dict` that maps variables into their default values
-    for the current-level system.
-  - `get_noiseeqs(sys)`: Noise equations of the current-level system.
-  - `get_description(sys)`: A string that describes what a system represents.
-  - `get_metadata(sys)`: Any metadata about the system or its origin to be used by downstream packages.
-
-Note that if you know a system is an `AbstractTimeDependentSystem` you could use `get_iv` to get the
-unique independent variable directly, rather than using `independent_variables(sys)[1]`, which is clunky and may cause problems if `sys` is an `AbstractMultivariateSystem` because there may be more than one independent variable. `AbstractTimeIndependentSystem`s do not have a method `get_iv`, and `independent_variables(sys)` will return a size-zero result for such. For an `AbstractMultivariateSystem`, `get_ivs` is equivalent.
-
-For the `parameters`, `unknowns`, `continuous_events`, and `discrete_events` accessors there are corresponding `parameters_toplevel`, `unknowns_toplevel`, `continuous_events_toplevel`, and `discrete_events_toplevel` accessors which work similarly, but ignore the content of subsystems. Furthermore, a `equations_toplevel` version of `equations` exists as well, however, it can only be applied to non-complete systems.
-
-A system could also have caches:
-
-  - `get_jac(sys)`: The Jacobian of a system.
-  - `get_tgrad(sys)`: The gradient with respect to time of a system.
-
-## Transformations
-
-Transformations are functions which send a valid `AbstractSystem` definition to
-another `AbstractSystem`. These are passes, like optimizations (e.g., Block-Lower
-Triangle transformations), or changes to the representation, which allow for
-alternative numerical methods to be utilized on the model (e.g., DAE index reduction).
-
-## Analyses
-
-Analyses are functions on a system which return information about the corresponding
-properties, like whether its parameters are structurally identifiable, or whether
-it's linear.
-
-## Function Calculation and Generation
-
-The calculation and generation functions allow for calculating additional
-quantities to enhance the numerical methods applied to the resulting system.
-The calculations, like `calculate_jacobian`, generate ModelingToolkit IR for
-the Jacobian of the system, while the generations, like `generate_jacobian`,
-generate compiled output for the numerical solvers by applying `build_function`
-to the generated code. Additionally, many systems have function-type outputs,
-which cobble together the generation functionality for a system, for example,
-`ODEFunction` can be used to generate a DifferentialEquations-based `ODEFunction`
-with compiled version of the ODE itself, the Jacobian, the mass matrix, etc.
-
-Below are the possible calculation and generation functions:
-
-```@docs
-calculate_tgrad
-calculate_gradient
-calculate_jacobian
-calculate_factorized_W
-calculate_hessian
-generate_tgrad
-generate_gradient
-generate_jacobian
-generate_factorized_W
-generate_hessian
-```
-
-Additionally, `jacobian_sparsity(sys)` and `hessian_sparsity(sys)`
-exist on the appropriate systems for fast generation of the sparsity
-patterns via an abstract interpretation without requiring differentiation.
-
-## Problem Constructors
-
-At the end, the system types have `DEProblem` constructors, like `ODEProblem`,
-which allow for directly generating the problem types required for numerical
-methods. The first argument is always the `AbstractSystem`, and the next
-arguments match the argument order of their original constructors. Whenever an
-array would normally be provided, such as `u0` the initial condition of an
-`ODEProblem`, it is instead replaced with a variable map, i.e., an array of
-pairs `var=>value`, which allows the user to designate the values without having
-to know the order that ModelingToolkit is internally using.
-
-For the value maps, the parameters are allowed to be functions of each other,
-and value maps of unknowns can be functions of the parameters, i.e. you can do:
-
-```
-u0 = [
-  lorenz1.x => 2.0
-  lorenz2.x => lorenz1.x * lorenz1.p
-]
-```
-
-## Default Value Handling
-
-The `AbstractSystem` types allow for specifying default values, for example
-`defaults` inside of them. At problem construction time, these values are merged
-into the value maps, where for any repeats the value maps override the default.
-In addition, defaults of a higher level in the system override the defaults of
-a lower level in the system.
-
-## Namespacing
-
-By default, unsimplified systems will namespace variables accessed via `getproperty`.
-Systems created via `@mtkcompile`, or ones passed through `mtkcompile` or
-`complete` will not perform this namespacing. However, all of these processes modify
-the system in a variety of ways. To toggle namespacing without transforming any other
-property of the system, use `toggle_namespacing`.
-
-```@docs
-toggle_namespacing
-```
diff --git a/docs/src/basics/Debugging.md b/docs/src/basics/Debugging.md
index 4f9c2c07d7..ccf1de263a 100644
--- a/docs/src/basics/Debugging.md
+++ b/docs/src/basics/Debugging.md
@@ -68,6 +68,6 @@ dprob[ModelingToolkit.ASSERTION_LOG_VARIABLE] = false;
 solve(dprob, Tsit5());
 ```
 
-```@docs
+```@docs; canonical = false
 debug_system
 ```
diff --git a/docs/src/basics/Events.md b/docs/src/basics/Events.md
index 89f874b08b..aca8eb4b68 100644
--- a/docs/src/basics/Events.md
+++ b/docs/src/basics/Events.md
@@ -249,11 +249,11 @@ par = @parameters g = 9.8
 bb_eqs = [D(x) ~ v
           D(v) ~ -g]
 
-function bb_affect!(integ, u, p, ctx)
-    integ.u[u.v] = -integ.u[u.v]
+function bb_affect!(mod, obs, integ, ctx)
+    return (; v = -mod.v)
 end
 
-reflect = [x ~ 0] => (bb_affect!, [v], [], [], nothing)
+reflect = [x ~ 0] => (bb_affect!, (; v))
 
 @mtkcompile bb_sys = System(bb_eqs, t, sts, par,
     continuous_events = reflect)
diff --git a/docs/src/basics/FAQ.md b/docs/src/basics/FAQ.md
index e3b12b46ab..1a1ffe75ca 100644
--- a/docs/src/basics/FAQ.md
+++ b/docs/src/basics/FAQ.md
@@ -88,7 +88,7 @@ Strings are not considered symbolic variables, and thus cannot directly be used
 indexing. However, ModelingToolkit does provide a method to parse the string representation of
 a variable, given the system in which that variable exists.
 
-```@docs
+```@docs; canonical = false
 ModelingToolkit.parse_variable
 ```
 
@@ -260,14 +260,14 @@ D = Differential(x)
 Tunable parameters are floating point parameters, not used in callbacks and not marked with `tunable = false` in their metadata. These are expected to be used with AD
 and optimization libraries. As such, they are stored together in one `Vector{T}`. To obtain the ordering of tunable parameters in this buffer, use:
 
-```@docs
+```@docs; canonical = false
 tunable_parameters
 ```
 
 If you have an array in which a particular dimension is in the order of tunable parameters (e.g. the jacobian with respect to tunables) then that dimension of the
 array can be reordered into the required permutation using the symbolic variables:
 
-```@docs
+```@docs; canonical = false
 reorder_dimension_by_tunables!
 reorder_dimension_by_tunables
 ```
diff --git a/docs/src/basics/InputOutput.md b/docs/src/basics/InputOutput.md
index 2e9da1c2db..b1eb2905df 100644
--- a/docs/src/basics/InputOutput.md
+++ b/docs/src/basics/InputOutput.md
@@ -75,7 +75,7 @@ u = [rand()]
 
 ## Generating an output function, ``g``
 
-ModelingToolkit can also generate a function that computes a specified output of a system, the function ``y = g(x, u, p, t)`` above. This is done using the function [`build_explicit_observed_function`](@ref). When generating an output function, the user must specify the output variable(s) of interest, as well as any inputs if inputs are relevant to compute the output.
+ModelingToolkit can also generate a function that computes a specified output of a system, the function ``y = g(x, u, p, t)`` above. This is done using the function [`ModelingToolkit.build_explicit_observed_function`](@ref). When generating an output function, the user must specify the output variable(s) of interest, as well as any inputs if inputs are relevant to compute the output.
 
 The order of the user-specified output variables determines the order of the output vector ``y``.
 
@@ -93,7 +93,7 @@ See [Linearization](@ref linearization).
 Pages = ["InputOutput.md"]
 ```
 
-```@docs
+```@docs; canonical=false
 ModelingToolkit.generate_control_function
 ModelingToolkit.build_explicit_observed_function
 ```
diff --git a/docs/src/basics/Linearization.md b/docs/src/basics/Linearization.md
index 3951aab28a..0d219d35a5 100644
--- a/docs/src/basics/Linearization.md
+++ b/docs/src/basics/Linearization.md
@@ -153,7 +153,7 @@ Also see [ControlSystemsMTK.jl](https://juliacontrol.github.io/ControlSystemsMTK
 Pages = ["Linearization.md"]
 ```
 
-```@docs
+```@docs; canonical = false
 linearize
 ModelingToolkit.linearize_symbolic
 ModelingToolkit.linearization_function
diff --git a/docs/src/internals.md b/docs/src/internals.md
index 0381e854b0..409ca51571 100644
--- a/docs/src/internals.md
+++ b/docs/src/internals.md
@@ -2,45 +2,3 @@
 
 This is a page for detailing some of the inner workings to help future
 contributors to the library.
-
-## Observables and Variable Elimination
-
-In the variable “elimination” algorithms, what is actually done is that variables
-are removed from being unknowns and equations are moved into the `observed` category
-of the system. The `observed` equations are explicit algebraic equations which
-are then substituted out to completely eliminate these variables from the other
-equations, allowing the system to act as though these variables no longer exist.
-
-However, a user may want to interact with such variables, for example,
-plotting their output. For this reason, these relationships are stored,
-and are then used to generate the `observed` equation found in the
-`SciMLFunction` interface, so that `sol[x]` lazily reconstructs the observed
-variable when necessary. In this sense, there is an equivalence between
-observables and the variable elimination system.
-
-The procedure for variable elimination inside [`mtkcompile`](@ref) is
-
- 1. [`ModelingToolkit.initialize_system_structure`](@ref).
- 2. [`ModelingToolkit.alias_elimination`](@ref). This step moves equations into `observed(sys)`.
- 3. [`ModelingToolkit.dae_index_lowering`](@ref) by means of [`pantelides!`](@ref) (if the system is an [`System`](@ref)).
- 4. [`ModelingToolkit.tearing`](@ref).
-
-## Preparing a system for simulation
-
-Before a simulation or optimization can be performed, the symbolic equations stored in an [`AbstractSystem`](@ref) must be converted into executable code. This step typically occurs after the simplification explained above, and is performed when an instance of a [`SciMLBase.AbstractSciMLProblem`](@ref), such as a [`ODEProblem`](@ref), is constructed.
-The call chain typically looks like this, with the function names in the case of an `System` indicated in parentheses
-
- 1. Problem constructor ([`ODEProblem`](@ref))
- 2. Build an `DEFunction` ([`process_DEProblem`](@ref) -> [`ODEFunction`](@ref)
- 3. Write actual executable code ([`generate_function`](@ref) or [`generate_custom_function`](@ref))
-
-Apart from [`generate_function`](@ref), which generates the dynamics function, `ODEFunction` also builds functions for observed equations (`build_explicit_observed_function`) and Jacobians (`generate_jacobian`) etc. These are all stored in the `ODEFunction`.
-
-## Creating an `MTKParameters` object
-
-It may be useful to create a parameter object without creating the problem. For this
-purpose, the `MTKParameters` constructor is exposed as public API.
-
-```@docs
-MTKParameters
-```
diff --git a/docs/src/tutorials/SampledData.md b/docs/src/tutorials/SampledData.md
index c2d6c9308d..9f3f340f46 100644
--- a/docs/src/tutorials/SampledData.md
+++ b/docs/src/tutorials/SampledData.md
@@ -189,7 +189,7 @@ connections = [r ~ sin(t)          # reference signal
 @named cl = System(connections, t, systems = [f, c, p])
 ```
 
-```@docs
+```@docs; canonical = false
 Sample
 Hold
 ShiftIndex
diff --git a/docs/src/tutorials/disturbance_modeling.md b/docs/src/tutorials/disturbance_modeling.md
index 341077b76b..cdbfa25b40 100644
--- a/docs/src/tutorials/disturbance_modeling.md
+++ b/docs/src/tutorials/disturbance_modeling.md
@@ -123,7 +123,7 @@ y &= g(x, u, p, t)
 \end{aligned}
 ```
 
-To make use of the model defined above for state estimation, we may want to generate a Julia function for the dynamics ``f`` and the output equations ``g`` that we can plug into, e.g., a nonlinear version of a Kalman filter or a particle filter, etc. MTK contains utilities to do this, namely, [`generate_control_function`](@ref) and [`build_explicit_observed_function`](@ref) (described in more details in ["Input output"](@ref inputoutput)). These functions take keyword arguments `disturbance_inputs` and `disturbance_argument`, that indicate which variables in the model are considered part of ``w``, and whether or not these variables are to be added as function arguments to ``f``, i.e., whether we have ``f(x, u, p, t)`` or ``f(x, u, p, t, w)``. If we do not include the disturbance inputs as function arguments, MTK will assume that the ``w`` variables are all zero, but any dynamics associated with these variables, such as disturbance models, will be included in the generated function. This allows a state estimator to estimate the state of the disturbance model, provided that this state is [observable](https://en.wikipedia.org/wiki/Observability) from the measured outputs of the system.
+To make use of the model defined above for state estimation, we may want to generate a Julia function for the dynamics ``f`` and the output equations ``g`` that we can plug into, e.g., a nonlinear version of a Kalman filter or a particle filter, etc. MTK contains utilities to do this, namely, [`ModelingToolkit.generate_control_function`](@ref) and [`ModelingToolkit.build_explicit_observed_function`](@ref) (described in more details in ["Input output"](@ref inputoutput)). These functions take keyword arguments `disturbance_inputs` and `disturbance_argument`, that indicate which variables in the model are considered part of ``w``, and whether or not these variables are to be added as function arguments to ``f``, i.e., whether we have ``f(x, u, p, t)`` or ``f(x, u, p, t, w)``. If we do not include the disturbance inputs as function arguments, MTK will assume that the ``w`` variables are all zero, but any dynamics associated with these variables, such as disturbance models, will be included in the generated function. This allows a state estimator to estimate the state of the disturbance model, provided that this state is [observable](https://en.wikipedia.org/wiki/Observability) from the measured outputs of the system.
 
 Below, we demonstrate
 
diff --git a/docs/src/tutorials/initialization.md b/docs/src/tutorials/initialization.md
index 90d86f1521..a804ca1b3e 100644
--- a/docs/src/tutorials/initialization.md
+++ b/docs/src/tutorials/initialization.md
@@ -362,7 +362,7 @@ which requires the system-level information and the additional nonlinear equatio
 tagged to the system.
 
 ```@example init
-isys = generate_initializesystem(pend, u0map = [x => 1.0, y => 0.0], guesses = [λ => 1])
+isys = generate_initializesystem(pend; op = [x => 1.0, y => 0.0], guesses = [λ => 1])
 ```
 
 We can inspect what its equations and unknown values are:
@@ -388,7 +388,7 @@ does not match the number of unknowns, which we can use to investigate our overd
 
 ```@example init
 isys = ModelingToolkit.generate_initializesystem(
-    pend, u0map = [x => 1, y => 0.0, D(y) => 2.0, λ => 1], guesses = [λ => 1])
+    pend; op = [x => 1, y => 0.0, D(y) => 2.0, λ => 1], guesses = [λ => 1])
 ```
 
 ```@example init
@@ -418,7 +418,7 @@ creates the special initialization system for a given `sys`. This is done as fol
 
 ```@example init
 iprob = ModelingToolkit.InitializationProblem(pend, 0.0,
-    [x => 1, y => 0.0, D(y) => 2.0, λ => 1], [g => 1], guesses = [λ => 1])
+    [x => 1, y => 0.0, D(y) => 2.0, λ => 1, g => 1], guesses = [λ => 1])
 ```
 
 We can see that because the system is overdetermined we receive a NonlinearLeastSquaresProblem,
@@ -460,7 +460,7 @@ some of the conditions:
 
 ```@example init
 iprob = ModelingToolkit.InitializationProblem(pend, 0.0,
-    [x => 1, y => 0.0, D(y) => 0.0, λ => 0], [g => 1], guesses = [λ => 1])
+    [x => 1, y => 0.0, D(y) => 0.0, λ => 0, g => 1], guesses = [λ => 1])
 ```
 
 gives a NonlinearLeastSquaresProblem which can be solved:
@@ -477,7 +477,7 @@ In comparison, if we have a well-conditioned system:
 
 ```@example init
 iprob = ModelingToolkit.InitializationProblem(pend, 0.0,
-    [x => 1, y => 0.0], [g => 1], guesses = [λ => 1])
+    [x => 1, y => 0.0, g => 1], guesses = [λ => 1])
 ```
 
 notice that we instead obtained a NonlinearSystem. In this case we have to use
diff --git a/docs/src/tutorials/linear_analysis.md b/docs/src/tutorials/linear_analysis.md
index 250dbaa6a7..2119ea9ad1 100644
--- a/docs/src/tutorials/linear_analysis.md
+++ b/docs/src/tutorials/linear_analysis.md
@@ -146,7 +146,7 @@ nyquistplot(P)
 Pages = ["linear_analysis.md"]
 ```
 
-```@autodocs
+```@autodocs; canonical = false
 Modules = [ModelingToolkit]
 Pages   = ["systems/analysis_points.jl"]
 Order   = [:function, :type]
diff --git a/docs/src/tutorials/stochastic_diffeq.md b/docs/src/tutorials/stochastic_diffeq.md
index 72a77eda05..d7326c36c6 100644
--- a/docs/src/tutorials/stochastic_diffeq.md
+++ b/docs/src/tutorials/stochastic_diffeq.md
@@ -39,7 +39,7 @@ By "multiplying" the equations by $dt$, the notation used in
 can be recovered.
 
 We use this Langevin-like notation because it allows us to extend MTK modeling capacity from ODEs to SDEs,
-using only a single new concept, `@brownian` variables, which represent $\frac{dB}{dt}$ in the above equation.
+using only a single new concept, `@brownians` variables, which represent $\frac{dB}{dt}$ in the above equation.
 
 ```@example SDE
 using ModelingToolkit, StochasticDiffEq
@@ -48,7 +48,7 @@ using Plots
 
 @parameters σ=10.0 ρ=2.33 β=26.0
 @variables x(t)=5.0 y(t)=5.0 z(t)=1.0
-@brownian B
+@brownians B
 eqs = [D(x) ~ σ * (y - x) + 0.3x * B,
     D(y) ~ x * (ρ - z) - y + 0.3y * B,
     D(z) ~ x * y - β * z + 0.3z * B]
@@ -80,10 +80,10 @@ plot(sol)
 ```
 
 If you want uncorrelated noise for each equation,
-multiple `@brownian` variables have to be declared.
+multiple `@brownians` variables have to be declared.
 
 ```@example SDE
-@brownian Bx By Bz
+@brownians Bx By Bz
 eqs = [D(x) ~ σ * (y - x) + 0.3x * Bx,
     D(y) ~ x * (ρ - z) - y + 0.3y * By,
     D(z) ~ x * y - β * z + 0.3z * Bz]
diff --git a/src/ModelingToolkit.jl b/src/ModelingToolkit.jl
index 255c5ade48..0a7fdd5209 100644
--- a/src/ModelingToolkit.jl
+++ b/src/ModelingToolkit.jl
@@ -330,12 +330,12 @@ export toexpr, get_variables
 export simplify, substitute
 export build_function
 export modelingtoolkitize
-export generate_initializesystem, Initial, isinitial
+export generate_initializesystem, Initial, isinitial, InitializationProblem
 
 export alg_equations, diff_equations, has_alg_equations, has_diff_equations
 export get_alg_eqs, get_diff_eqs, has_alg_eqs, has_diff_eqs
 
-export @variables, @parameters, @independent_variables, @constants, @brownian
+export @variables, @parameters, @independent_variables, @constants, @brownians, @brownian
 export @named, @nonamespace, @namespace, extend, compose, complete, toggle_namespacing
 export debug_system
 
@@ -361,6 +361,16 @@ export AbstractCollocation, JuMPCollocation, InfiniteOptCollocation,
        CasADiCollocation, PyomoCollocation
 export DynamicOptSolution
 
-@public apply_to_variables
+@public apply_to_variables, equations_toplevel, unknowns_toplevel, parameters_toplevel
+@public continuous_events_toplevel, discrete_events_toplevel, assertions, is_alg_equation
+@public is_diff_equation, Equality, linearize_symbolic, reorder_unknowns
+@public similarity_transform, inputs, outputs, bound_inputs, unbound_inputs, bound_outputs
+@public unbound_outputs, is_bound
+
+for prop in [SYS_PROPS; [:continuous_events, :discrete_events]]
+    getter = Symbol(:get_, prop)
+    hasfn = Symbol(:has_, prop)
+    @eval @public $getter, $hasfn
+end
 
 end # module
diff --git a/src/deprecations.jl b/src/deprecations.jl
index aaafa0d283..56c166f575 100644
--- a/src/deprecations.jl
+++ b/src/deprecations.jl
@@ -170,3 +170,10 @@ for T in [:NonlinearProblem, :NonlinearLeastSquaresProblem,
         end
     end
 end
+
+macro brownian(xs...)
+    return quote
+        Base.depwarn("`@brownian` is deprecated. Use `@brownians` instead", :brownian_macro)
+        $(@__MODULE__).@brownians $(xs...)
+    end
+end
diff --git a/src/discretedomain.jl b/src/discretedomain.jl
index 73e6b8b7fa..a10fec81b4 100644
--- a/src/discretedomain.jl
+++ b/src/discretedomain.jl
@@ -1,5 +1,11 @@
 using Symbolics: Operator, Num, Term, value, recursive_hasoperator
 
+"""
+    function SampleTime()
+
+`SampleTime()` can be used in the equations of a hybrid system to represent time sampled
+at the inferred clock for that equation.
+"""
 struct SampleTime <: Operator
     SampleTime() = SymbolicUtils.term(SampleTime, type = Real)
 end
diff --git a/src/problems/compatibility.jl b/src/problems/compatibility.jl
index 36302a9e74..e5608ce4f1 100644
--- a/src/problems/compatibility.jl
+++ b/src/problems/compatibility.jl
@@ -119,7 +119,7 @@ function check_has_noise(sys::System, T)
         """
         if !isempty(brownians(sys))
             msg = """
-            Systems constructed by defining Brownian variables with `@brownian` must be \
+            Systems constructed by defining Brownian variables with `@brownians` must be \
             simplified by calling `mtkcompile` before a `$T` can be constructed.
             """
         end
diff --git a/src/structural_transformation/symbolics_tearing.jl b/src/structural_transformation/symbolics_tearing.jl
index af0973f1dc..f8f05ffb7f 100644
--- a/src/structural_transformation/symbolics_tearing.jl
+++ b/src/structural_transformation/symbolics_tearing.jl
@@ -1265,7 +1265,7 @@ function tearing_hacks(sys, obs, unknowns, neweqs; array = true)
     for (arrvar, cnt) in arr_obs_occurrences
         cnt == length(arrvar) || continue
         # firstindex returns 1 for multidimensional array symbolics
-        firstind = first(eachindex(arrvar))
+        firstind = Tuple(first(eachindex(arrvar)))
         scal = [arrvar[i] for i in eachindex(arrvar)]
         # respect non-1-indexed arrays
         # TODO: get rid of this hack together with the above hack, then remove OffsetArrays dependency
@@ -1273,8 +1273,7 @@ function tearing_hacks(sys, obs, unknowns, neweqs; array = true)
         # try to `create_array(OffsetArray{...}, ...)` which errors.
         # `term(Origin(firstind), scal)` doesn't retain the `symtype` and `size`
         # of `scal`.
-        rhs = scal
-        rhs = change_origin(firstind, rhs)
+        rhs = change_origin(firstind, scal)
         push!(obs_arr_eqs, arrvar ~ rhs)
     end
     append!(obs, obs_arr_eqs)
@@ -1284,7 +1283,7 @@ end
 
 # PART OF HACK
 function change_origin(origin, arr)
-    if all(isone, Tuple(origin))
+    if all(isone, origin)
         return arr
     end
     return Origin(origin)(arr)
diff --git a/src/systems/abstractsystem.jl b/src/systems/abstractsystem.jl
index 941739f127..386dd2f2fd 100644
--- a/src/systems/abstractsystem.jl
+++ b/src/systems/abstractsystem.jl
@@ -730,47 +730,49 @@ function unflatten_parameters!(buffer, params, all_ps)
     end
 end
 
-for prop in [:eqs
-             :tag
-             :noise_eqs
-             :iv
-             :unknowns
-             :ps
-             :tspan
-             :brownians
-             :jumps
-             :name
-             :description
-             :var_to_name
-             :defaults
-             :guesses
-             :observed
-             :systems
-             :constraints
-             :bcs
-             :domain
-             :ivs
-             :dvs
-             :connector_type
-             :preface
-             :initializesystem
-             :initialization_eqs
-             :schedule
-             :tearing_state
-             :metadata
-             :gui_metadata
-             :is_initializesystem
-             :is_discrete
-             :parameter_dependencies
-             :assertions
-             :ignored_connections
-             :parent
-             :is_dde
-             :tstops
-             :index_cache
-             :isscheduled
-             :costs
-             :consolidate]
+const SYS_PROPS = [:eqs
+                   :tag
+                   :noise_eqs
+                   :iv
+                   :unknowns
+                   :ps
+                   :tspan
+                   :brownians
+                   :jumps
+                   :name
+                   :description
+                   :var_to_name
+                   :defaults
+                   :guesses
+                   :observed
+                   :systems
+                   :constraints
+                   :bcs
+                   :domain
+                   :ivs
+                   :dvs
+                   :connector_type
+                   :preface
+                   :initializesystem
+                   :initialization_eqs
+                   :schedule
+                   :tearing_state
+                   :metadata
+                   :gui_metadata
+                   :is_initializesystem
+                   :is_discrete
+                   :parameter_dependencies
+                   :assertions
+                   :ignored_connections
+                   :parent
+                   :is_dde
+                   :tstops
+                   :index_cache
+                   :isscheduled
+                   :costs
+                   :consolidate]
+
+for prop in SYS_PROPS
     fname_get = Symbol(:get_, prop)
     fname_has = Symbol(:has_, prop)
     @eval begin
@@ -780,7 +782,6 @@ for prop in [:eqs
         Get the internal field `$($(QuoteNode(prop)))` of a system `sys`.
         It only includes `$($(QuoteNode(prop)))` local to `sys`; not those of its subsystems,
         like `unknowns(sys)`, `parameters(sys)` and `equations(sys)` does.
-        This is equivalent to, but preferred over `sys.$($(QuoteNode(prop)))`.
 
         See also [`has_$($(QuoteNode(prop)))`](@ref).
         """
@@ -1519,7 +1520,7 @@ $(TYPEDSIGNATURES)
 Get the default values of the system sys and its subsystems.
 If they are not explicitly provided, variables and parameters are initialized to these values.
 
-See also [`initialization_equations`](@ref), [`parameter_dependencies`](@ref) and [`ModelingToolkit.get_defaults`](@ref).
+See also [`initialization_equations`](@ref) and [`ModelingToolkit.get_defaults`](@ref).
 """
 function defaults(sys::AbstractSystem)
     systems = get_systems(sys)
@@ -1722,7 +1723,7 @@ $(TYPEDSIGNATURES)
 
 Get the initialization equations of the system `sys` and its subsystems.
 
-See also [`guesses`](@ref), [`defaults`](@ref), [`parameter_dependencies`](@ref) and [`ModelingToolkit.get_initialization_eqs`](@ref).
+See also [`guesses`](@ref), [`defaults`](@ref) and [`ModelingToolkit.get_initialization_eqs`](@ref).
 """
 function initialization_equations(sys::AbstractSystem)
     eqs = get_initialization_eqs(sys)
diff --git a/src/systems/callbacks.jl b/src/systems/callbacks.jl
index a12dad9d72..cf9c53e610 100644
--- a/src/systems/callbacks.jl
+++ b/src/systems/callbacks.jl
@@ -941,7 +941,23 @@ function discrete_events(sys::AbstractSystem)
     cbs
 end
 
+"""
+    $(TYPEDSIGNATURES)
+
+Returns whether the system `sys` has the internal field `discrete_events`.
+
+See also [`get_discrete_events`](@ref).
+"""
 has_discrete_events(sys::AbstractSystem) = isdefined(sys, :discrete_events)
+"""
+    $(TYPEDSIGNATURES)
+
+Get the internal field `discrete_events` of a system `sys`.
+It only includes `discrete_events` local to `sys`; not those of its subsystems,
+like `unknowns(sys)`, `parameters(sys)` and `equations(sys)` does.
+
+See also [`has_discrete_events`](@ref).
+"""
 function get_discrete_events(sys::AbstractSystem)
     has_discrete_events(sys) || return SymbolicDiscreteCallback[]
     getfield(sys, :discrete_events)
@@ -984,7 +1000,23 @@ function continuous_events(sys::AbstractSystem)
     filter(!isempty, cbs)
 end
 
+"""
+    $(TYPEDSIGNATURES)
+
+Returns whether the system `sys` has the internal field `continuous_events`.
+
+See also [`get_continuous_events`](@ref).
+"""
 has_continuous_events(sys::AbstractSystem) = isdefined(sys, :continuous_events)
+"""
+    $(TYPEDSIGNATURES)
+
+Get the internal field `continuous_events` of a system `sys`.
+It only includes `continuous_events` local to `sys`; not those of its subsystems,
+like `unknowns(sys)`, `parameters(sys)` and `equations(sys)` does.
+
+See also [`has_continuous_events`](@ref).
+"""
 function get_continuous_events(sys::AbstractSystem)
     has_continuous_events(sys) || return SymbolicContinuousCallback[]
     getfield(sys, :continuous_events)
diff --git a/src/systems/codegen.jl b/src/systems/codegen.jl
index 17ff652c41..96ff14f58a 100644
--- a/src/systems/codegen.jl
+++ b/src/systems/codegen.jl
@@ -277,6 +277,16 @@ function generate_tgrad(
         expression, wrap_gfw, (2, 3, is_split(sys)), res; eval_expression, eval_module)
 end
 
+"""
+    $(TYPEDSIGNATURES)
+
+Return an array of symbolic hessians corresponding to the equations of the system.
+
+# Keyword Arguments
+
+- `sparse`: Controls whether the symbolic hessians are sparse matrices
+- `simplify`: Forwarded to `Symbolics.hessian`
+"""
 function calculate_hessian(sys::System; simplify = false, sparse = false)
     rhs = [eq.rhs - eq.lhs for eq in full_equations(sys)]
     dvs = unknowns(sys)
@@ -484,6 +494,17 @@ function W_sparsity(sys::System)
     jac_sparsity .| M_sparsity
 end
 
+"""
+    $(TYPEDSIGNATURES)
+
+Return the matrix to use as the jacobian prototype given the W-sparsity matrix of the
+system. This is not the same as the jacobian sparsity pattern.
+
+# Keyword arguments
+
+- `u0`: The `u0` vector for the problem.
+- `sparse`: The prototype is `nothing` for non-sparse matrices.
+"""
 function calculate_W_prototype(W_sparsity; u0 = nothing, sparse = false)
     sparse || return nothing
     uElType = u0 === nothing ? Float64 : eltype(u0)
@@ -599,6 +620,12 @@ function generate_cost(sys::System; expression = Val{true}, wrap_gfw = Val{false
         expression, wrap_gfw, (2, nargs, is_split(sys)), res; eval_expression, eval_module)
 end
 
+"""
+    $(TYPEDSIGNATURES)
+
+Calculate the gradient of the consolidated cost of `sys` with respect to the unknowns.
+`simplify` is forwarded to `Symbolics.gradient`.
+"""
 function calculate_cost_gradient(sys::System; simplify = false)
     obj = cost(sys)
     dvs = unknowns(sys)
@@ -629,6 +656,13 @@ function generate_cost_gradient(
         expression, wrap_gfw, (2, 2, is_split(sys)), res; eval_expression, eval_module)
 end
 
+"""
+    $(TYPEDSIGNATURES)
+
+Calculate the hessian of the consolidated cost of `sys` with respect to the unknowns.
+`simplify` is forwarded to `Symbolics.hessian`. `sparse` controls whether a sparse
+matrix is returned.
+"""
 function calculate_cost_hessian(sys::System; sparse = false, simplify = false)
     obj = cost(sys)
     dvs = unknowns(sys)
diff --git a/src/systems/connectors.jl b/src/systems/connectors.jl
index cf2f24a7d1..5d6227d4c1 100644
--- a/src/systems/connectors.jl
+++ b/src/systems/connectors.jl
@@ -860,6 +860,12 @@ function expand_variable_connections(sys::AbstractSystem; ignored_variables = no
     return sys
 end
 
+"""
+    function expand_connections(sys::AbstractSystem)
+
+Given a hierarchical system with [`connect`](@ref) equations, expand the connection
+equations and return the new system.
+"""
 function expand_connections(sys::AbstractSystem;
         debug = false, tol = 1e-10, scalarize = true)
     sys = remove_analysis_points(sys)
diff --git a/src/systems/diffeqs/basic_transformations.jl b/src/systems/diffeqs/basic_transformations.jl
index 430260f60a..2eeb1a05d5 100644
--- a/src/systems/diffeqs/basic_transformations.jl
+++ b/src/systems/diffeqs/basic_transformations.jl
@@ -483,7 +483,7 @@ function noise_to_brownians(sys::System; names::Union{Symbol, Vector{Symbol}} =
         """))
     end
     brownvars = map(names) do name
-        unwrap(only(@brownian $name))
+        unwrap(only(@brownians $name))
     end
 
     terms = if ndims(neqs) == 1
diff --git a/src/systems/nonlinear/initializesystem.jl b/src/systems/nonlinear/initializesystem.jl
index 5dfd69b72f..4f7bc7771f 100644
--- a/src/systems/nonlinear/initializesystem.jl
+++ b/src/systems/nonlinear/initializesystem.jl
@@ -1,3 +1,34 @@
+"""
+    $(TYPEDSIGNATURES)
+
+Generate the initialization system for `sys`. The initialization system is a system of
+nonlinear equations that solve for the full set of initial conditions of `sys` given
+specified constraints.
+
+The initialization system can be of two types: time-dependent and time-independent.
+Time-dependent initialization systems solve for the initial values of unknowns as well as
+the values of solvable parameters of the system. Time-independent initialization systems
+only solve for solvable parameters of the system.
+
+# Keyword arguments
+
+- `time_dependent_init`: Whether to create an initialization system for a time-dependent
+  system. A time-dependent initialization requires a time-dependent `sys`, but a time-
+  independent initialization can be created regardless.
+- `op`: The operating point of user-specified initial conditions of variables in `sys`.
+- `initialization_eqs`: Additional initialization equations to use apart from those in
+  `initialization_equations(sys)`.
+- `guesses`: Additional guesses to use apart from those in `guesses(sys)`.
+- `default_dd_guess`: Default guess for dummy derivative variables in time-dependent
+  initialization.
+- `algebraic_only`: If `false`, does not use initialization equations (provided via the
+  keyword or part of the system) to construct initialization.
+- `check_defguess`: Whether to error when a variable does not have a default or guess
+  despite ModelingToolkit expecting it to.
+- `name`: The name of the initialization system.
+
+All other keyword arguments are forwarded to the [`System`](@ref) constructor.
+"""
 function generate_initializesystem(
         sys::AbstractSystem; time_dependent_init = is_time_dependent(sys), kwargs...)
     if time_dependent_init
diff --git a/src/variables.jl b/src/variables.jl
index 104616f6e4..38faf6dbda 100644
--- a/src/variables.jl
+++ b/src/variables.jl
@@ -163,18 +163,50 @@ function isvarkind(m, x)
     getmetadata(x, m, false)
 end
 
+"""
+    $(TYPEDSIGNATURES)
+
+Set the `input` metadata of variable `x` to `v`.
+"""
 setinput(x, v::Bool) = setmetadata(x, VariableInput, v)
+"""
+    $(TYPEDSIGNATURES)
+
+Set the `output` metadata of variable `x` to `v`.
+"""
 setoutput(x, v::Bool) = setmetadata(x, VariableOutput, v)
 setio(x, i::Bool, o::Bool) = setoutput(setinput(x, i), o)
 
+"""
+    $(TYPEDSIGNATURES)
+
+Check if variable `x` is marked as an input.
+"""
 isinput(x) = isvarkind(VariableInput, x)
+"""
+    $(TYPEDSIGNATURES)
+
+Check if variable `x` is marked as an output.
+"""
 isoutput(x) = isvarkind(VariableOutput, x)
 
 # Before the solvability check, we already have handled IO variables, so
 # irreducibility is independent from IO.
+"""
+    $(TYPEDSIGNATURES)
+
+Check if `x` is marked as irreducible. This prevents it from being eliminated as an
+observed variable in `mtkcompile`.
+"""
 isirreducible(x) = isvarkind(VariableIrreducible, x)
 setirreducible(x, v::Bool) = setmetadata(x, VariableIrreducible, v)
 state_priority(x::Union{Num, Symbolics.Arr}) = state_priority(unwrap(x))
+"""
+    $(TYPEDSIGNATURES)
+
+Return the `state_priority` metadata of variable `x`. This influences its priority to be
+chosen as a state in `mtkcompile`.
+"""
 state_priority(x) = convert(Float64, getmetadata(x, VariableStatePriority, 0.0))::Float64
 
 normalize_to_differential(x) = x
@@ -419,6 +451,11 @@ function getdescription(x)
     Symbolics.getmetadata(x, VariableDescription, "")
 end
 
+"""
+    $(TYPEDSIGNATURES)
+
+Check if variable `x` has a non-empty attached description.
+"""
 function hasdescription(x)
     getdescription(x) != ""
 end
@@ -438,11 +475,11 @@ $(SIGNATURES)
 
 Define one or more Brownian variables.
 """
-macro brownian(xs...)
+macro brownians(xs...)
     all(
         x -> x isa Symbol || Meta.isexpr(x, :call) && x.args[1] == :$ || Meta.isexpr(x, :$),
         xs) ||
-        error("@brownian only takes scalar expressions!")
+        error("@brownians only takes scalar expressions!")
     Symbolics._parse_vars(:brownian,
         Real,
         xs,
diff --git a/test/dde.jl b/test/dde.jl
index fb92bf61eb..df4acf377a 100644
--- a/test/dde.jl
+++ b/test/dde.jl
@@ -78,7 +78,7 @@ sol = solve(prob, RKMil(), seed = 100)
 
 @variables x(..) delx(t)
 @parameters a=-4.0 b=-2.0 c=10.0 α=-1.3 β=-1.2 γ=1.1
-@brownian η
+@brownians η
 τ = 1.0
 eqs = [D(x(t)) ~ a * x(t) + b * x(t - τ) + c + (α * x(t) + γ) * η, delx ~ x(t - τ)]
 @mtkcompile sys = System(eqs, t)
@@ -188,7 +188,7 @@ end
     alg = MethodOfSteps(Vern7())
     @test_nowarn solve(prob, alg)
 
-    @brownian r
+    @brownians r
     eqs = [D(x(t)) ~ -w * x(t - τ) + r]
     @named sys = System(eqs, t)
     sys = mtkcompile(sys)
diff --git a/test/debugging.jl b/test/debugging.jl
index 0cf80fd494..0ff7a0fb4d 100644
--- a/test/debugging.jl
+++ b/test/debugging.jl
@@ -3,7 +3,7 @@ import Logging
 using ModelingToolkit: t_nounits as t, D_nounits as D, ASSERTION_LOG_VARIABLE
 
 @variables x(t)
-@brownian a
+@brownians a
 @named inner_ode = System(D(x) ~ -sqrt(x), t; assertions = [(x > 0) => "ohno"])
 @named inner_sde = System([D(x) ~ -10sqrt(x) + 0.01a], t; assertions = [(x > 0) => "ohno"])
 sys_ode = mtkcompile(inner_ode)
diff --git a/test/initializationsystem.jl b/test/initializationsystem.jl
index 99e73e5b17..e12d05479f 100644
--- a/test/initializationsystem.jl
+++ b/test/initializationsystem.jl
@@ -592,7 +592,7 @@ end
 @testset "Initialization of parameters" begin
     @variables _x(..) y(t)
     @parameters p q
-    @brownian a b
+    @brownians a b
     x = _x(t)
     sarray_ctor = splat(SVector)
     # `System` constructor creates appropriate type with mtkcompile
@@ -881,7 +881,7 @@ end
 @testset "Update initializeprob parameters" begin
     @variables _x(..) y(t)
     @parameters p q
-    @brownian a b
+    @brownians a b
     x = _x(t)
 
     @testset "$Problem with $(SciMLBase.parameterless_type(typeof(alg)))" for (System, Problem, alg, rhss) in [
@@ -910,7 +910,7 @@ end
 @testset "Equations for dependent parameters" begin
     @variables _x(..)
     @parameters p q=5 r
-    @brownian a
+    @brownians a
     x = _x(t)
 
     @testset "$Problem with $(SciMLBase.parameterless_type(typeof(alg)))" for (System, Problem, alg, rhss) in [
@@ -933,7 +933,7 @@ end
 @testset "Re-creating initialization problem on remake" begin
     @variables _x(..) y(t)
     @parameters p q
-    @brownian a b
+    @brownians a b
     x = _x(t)
 
     @testset "$Problem with $(SciMLBase.parameterless_type(typeof(alg)))" for (Problem, alg, rhss) in [
@@ -965,7 +965,7 @@ end
 @testset "`remake` changes initialization problem types" begin
     @variables _x(..) y(t) z(t)
     @parameters p q
-    @brownian a
+    @brownians a
     x = _x(t)
 
     @testset "$Problem with $(SciMLBase.parameterless_type(typeof(alg)))" for (System, Problem, alg, rhss) in [
@@ -1022,7 +1022,7 @@ end
 @testset "`remake` preserves old u0map and pmap" begin
     @variables _x(..) y(t)
     @parameters p
-    @brownian a
+    @brownians a
     x = _x(t)
 
     @testset "$Problem with $(SciMLBase.parameterless_type(typeof(alg)))" for (System, Problem, alg, rhss) in [
@@ -1426,7 +1426,7 @@ end
 
 @testset "Trivial initialization is run on problem construction" begin
     @variables _x(..) y(t)
-    @brownian a
+    @brownians a
     @parameters tot
     x = _x(t)
     @testset "$Problem" for (Problem, lhs, rhs) in [
diff --git a/test/sciml_problem_inputs.jl b/test/sciml_problem_inputs.jl
index cd081118fb..46dfbbed21 100644
--- a/test/sciml_problem_inputs.jl
+++ b/test/sciml_problem_inputs.jl
@@ -150,7 +150,7 @@ end
 @testset "Deprecations" begin
     @variables _x(..) = 1.0
     @parameters p = 1.0
-    @brownian a
+    @brownians a
     x = _x(t)
     k = ShiftIndex(t)
 
diff --git a/test/sdesystem.jl b/test/sdesystem.jl
index 20f62c402c..f60e3a9777 100644
--- a/test/sdesystem.jl
+++ b/test/sdesystem.jl
@@ -591,7 +591,7 @@ end
 
 sts = @variables x(tt) y(tt) z(tt)
 ps = @parameters σ ρ
-@brownian β η
+@brownians β η
 s = 0.001
 β *= s
 η *= s
@@ -636,7 +636,7 @@ ssys = SDESystem(eqs, noise_eqs, t, [X], [p, d]; name = :ssys)
 # Issue#2814
 @parameters p d
 @variables x(tt)
-@brownian a
+@brownians a
 eqs = [D(x) ~ p - d * x + a * sqrt(p)]
 @mtkcompile sys = System(eqs, tt)
 u0 = @SVector[x => 10.0]
@@ -649,7 +649,7 @@ sprob = SDEProblem(sys, [u0; ps], tspan)
 
 # Ensure diagonal noise generates vector noise function
 @variables y(tt)
-@brownian b
+@brownians b
 eqs = [D(x) ~ p - d * x + a * sqrt(p)
        D(y) ~ p - d * y + b * sqrt(d)]
 @mtkcompile sys = System(eqs, tt)
@@ -663,7 +663,7 @@ sprob = SDEProblem(sys, [u0; ps], tspan)
 let
     @parameters σ ρ β
     @variables x(t) y(t) z(t)
-    @brownian a
+    @brownians a
     eqs = [D(x) ~ σ * (y - x) + 0.1a * x,
         D(y) ~ x * (ρ - z) - y + 0.1a * y,
         D(z) ~ x * y - β * z + 0.1a * z]
@@ -688,7 +688,7 @@ end
 
 let # test to make sure that scalar noise always receive the same kicks
     @variables x(t) y(t)
-    @brownian a
+    @brownians a
     eqs = [D(x) ~ a,
         D(y) ~ a]
 
@@ -701,7 +701,7 @@ end
 let # test that diagonal noise is correctly handled
     @parameters σ ρ β
     @variables x(t) y(t) z(t)
-    @brownian a b c
+    @brownians a b c
     eqs = [D(x) ~ σ * (y - x) + 0.1a * x,
         D(y) ~ x * (ρ - z) - y + 0.1b * y,
         D(z) ~ x * y - β * z + 0.1c * z]
@@ -728,7 +728,7 @@ end
 @testset "Non-diagonal noise check" begin
     @parameters σ ρ β
     @variables x(tt) y(tt) z(tt)
-    @brownian a b c d e f
+    @brownians a b c d e f
     eqs = [D(x) ~ σ * (y - x) + 0.1a * x + d,
         D(y) ~ x * (ρ - z) - y + 0.1b * y + e,
         D(z) ~ x * y - β * z + 0.1c * z + f]
@@ -757,7 +757,7 @@ end
 @testset "Diagonal noise, less brownians than equations" begin
     @parameters σ ρ β
     @variables x(tt) y(tt) z(tt)
-    @brownian a b
+    @brownians a b
     eqs = [D(x) ~ σ * (y - x) + 0.1a * x,     # One brownian
         D(y) ~ x * (ρ - z) - y + 0.1b * y, # Another brownian
         D(z) ~ x * y - β * z]              # no brownians -- still diagonal
@@ -781,7 +781,7 @@ end
 
 @testset "Passing `nothing` to `u0`" begin
     @variables x(t) = 1
-    @brownian b
+    @brownians b
     @mtkcompile sys = System([D(x) ~ x + b], t)
     prob = @test_nowarn SDEProblem(sys, nothing, (0.0, 1.0))
     @test_nowarn solve(prob, ImplicitEM())
@@ -794,7 +794,7 @@ end
         [input = true]
     end
     ps = @parameters a = 2
-    browns = @brownian η
+    browns = @brownians η
 
     eqs = [D(x) ~ -a * x + (input + 1) * η
            input ~ 0.0]
@@ -808,7 +808,7 @@ end
 
 @testset "Observed variables retained after `mtkcompile`" begin
     @variables x(t) y(t) z(t)
-    @brownian a
+    @brownians a
     @mtkcompile sys = System([D(x) ~ x + a, D(y) ~ y + a, z ~ x + y], t)
     @test length(observed(sys)) == 1
     prob = SDEProblem(sys, [x => 1.0, y => 1.0], (0.0, 1.0))
@@ -875,7 +875,7 @@ end
 @testset "Validate input types" begin
     @parameters p d
     @variables X(t)::Int64
-    @brownian z
+    @brownians z
     eq2 = D(X) ~ p - d * X + z
     @test_throws ModelingToolkit.ContinuousOperatorDiscreteArgumentError @mtkcompile ssys = System(
         [eq2], t)
@@ -929,7 +929,7 @@ end
 @testset "Error when constructing SDEProblem without `mtkcompile`" begin
     @parameters σ ρ β
     @variables x(tt) y(tt) z(tt)
-    @brownian a
+    @brownians a
     eqs = [D(x) ~ σ * (y - x) + 0.1a * x,
         D(y) ~ x * (ρ - z) - y + 0.1a * y,
         D(z) ~ x * y - β * z + 0.1a * z]
@@ -945,3 +945,7 @@ end
     de = mtkcompile(de)
     @test SDEProblem(de, [u0map; parammap], (0.0, 100.0)) isa SDEProblem
 end
+
+@testset "`@brownian` is deprecated" begin
+    @test_deprecated @brownian a b c
+end
diff --git a/test/structural_transformation/utils.jl b/test/structural_transformation/utils.jl
index 4c1ba1e807..4a2df411a6 100644
--- a/test/structural_transformation/utils.jl
+++ b/test/structural_transformation/utils.jl
@@ -358,7 +358,7 @@ end
     @testset "SDESystem" begin
         @variables x(t) p a
         @parameters y(t) q b
-        @brownian c
+        @brownians c
         @mtkcompile sys = System([D(x) ~ x + q * a, D(y) ~ y + p * b + c], t, [x, y],
             [p, q], [a, b, c]; initialization_eqs = [p + q ~ 4],
             guesses = [p => 1.0], defaults = [p => missing])
diff --git a/test/test_variable_metadata.jl b/test/test_variable_metadata.jl
index 482ebf927c..596372d45b 100644
--- a/test/test_variable_metadata.jl
+++ b/test/test_variable_metadata.jl
@@ -224,5 +224,5 @@ x = ModelingToolkit.toparam(x)
 @parameters y
 @test ModelingToolkit.getvariabletype(y) == ModelingToolkit.PARAMETER
 
-@brownian z
+@brownians z
 @test ModelingToolkit.getvariabletype(z) == ModelingToolkit.BROWNIAN