diff --git a/.gitignore b/.gitignore
index 3401a5a4b3..6d305ad348 100644
--- a/.gitignore
+++ b/.gitignore
@@ -4,3 +4,5 @@
 Manifest.toml
 .vscode
 .vscode/*
+docs/src/assets/Project.toml
+docs/src/assets/Manifest.toml
diff --git a/Project.toml b/Project.toml
index 470adc92b6..8cceafe504 100644
--- a/Project.toml
+++ b/Project.toml
@@ -1,7 +1,7 @@
 name = "ModelingToolkit"
 uuid = "961ee093-0014-501f-94e3-6117800e7a78"
 authors = ["Yingbo Ma <mayingbo5@gmail.com>", "Chris Rackauckas <accounts@chrisrackauckas.com> and contributors"]
-version = "8.71.1"
+version = "8.72.2"
 
 [deps]
 AbstractTrees = "1520ce14-60c1-5f80-bbc7-55ef81b5835c"
@@ -51,9 +51,11 @@ UnPack = "3a884ed6-31ef-47d7-9d2a-63182c4928ed"
 Unitful = "1986cc42-f94f-5a68-af5c-568840ba703d"
 
 [weakdeps]
+BifurcationKit = "0f109fa4-8a5d-4b75-95aa-f515264e7665"
 DeepDiffs = "ab62b9b5-e342-54a8-a765-a90f495de1a6"
 
 [extensions]
+MTKBifurcationKitExt = "BifurcationKit"
 MTKDeepDiffsExt = "DeepDiffs"
 
 [compat]
@@ -100,6 +102,7 @@ julia = "1.6"
 
 [extras]
 AmplNLWriter = "7c4d4715-977e-5154-bfe0-e096adeac482"
+BifurcationKit = "0f109fa4-8a5d-4b75-95aa-f515264e7665"
 BenchmarkTools = "6e4b80f9-dd63-53aa-95a3-0cdb28fa8baf"
 ControlSystemsMTK = "687d7614-c7e5-45fc-bfc3-9ee385575c88"
 DeepDiffs = "ab62b9b5-e342-54a8-a765-a90f495de1a6"
@@ -123,4 +126,4 @@ Sundials = "c3572dad-4567-51f8-b174-8c6c989267f4"
 Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
 
 [targets]
-test = ["AmplNLWriter", "BenchmarkTools", "ControlSystemsMTK", "NonlinearSolve", "ForwardDiff", "Ipopt", "Ipopt_jll", "ModelingToolkitStandardLibrary", "Optimization", "OptimizationOptimJL", "OptimizationMOI", "Random", "ReferenceTests", "SafeTestsets", "StableRNGs", "Statistics", "SteadyStateDiffEq", "Test", "StochasticDiffEq", "Sundials", "StochasticDelayDiffEq"]
+test = ["AmplNLWriter", "BifurcationKit", "BenchmarkTools", "ControlSystemsMTK", "NonlinearSolve", "ForwardDiff", "Ipopt", "Ipopt_jll", "ModelingToolkitStandardLibrary", "Optimization", "OptimizationOptimJL", "OptimizationMOI", "Random", "ReferenceTests", "SafeTestsets", "StableRNGs", "Statistics", "SteadyStateDiffEq", "Test", "StochasticDiffEq", "Sundials", "StochasticDelayDiffEq"]
diff --git a/docs/Project.toml b/docs/Project.toml
index 81a99e004e..68d33cc4c5 100644
--- a/docs/Project.toml
+++ b/docs/Project.toml
@@ -1,5 +1,6 @@
 [deps]
 BenchmarkTools = "6e4b80f9-dd63-53aa-95a3-0cdb28fa8baf"
+BifurcationKit = "0f109fa4-8a5d-4b75-95aa-f515264e7665"
 DifferentialEquations = "0c46a032-eb83-5123-abaf-570d42b7fbaa"
 Distributions = "31c24e10-a181-5473-b8eb-7969acd0382f"
 Documenter = "e30172f5-a6a5-5a46-863b-614d45cd2de4"
diff --git a/docs/pages.jl b/docs/pages.jl
index 0b6d0555cd..9d49c477ec 100644
--- a/docs/pages.jl
+++ b/docs/pages.jl
@@ -5,8 +5,10 @@ pages = [
         "tutorials/nonlinear.md",
         "tutorials/optimization.md",
         "tutorials/modelingtoolkitize.md",
+        "tutorials/programmatically_generating.md",
         "tutorials/stochastic_diffeq.md",
         "tutorials/parameter_identifiability.md",
+        "tutorials/bifurcation_diagram_computation.md",
         "tutorials/domain_connections.md"],
     "Examples" => Any["Basic Examples" => Any["examples/higher_order.md",
             "examples/spring_mass.md",
diff --git a/docs/src/basics/ContextualVariables.md b/docs/src/basics/ContextualVariables.md
index dcaebe95d7..de2802ce40 100644
--- a/docs/src/basics/ContextualVariables.md
+++ b/docs/src/basics/ContextualVariables.md
@@ -48,7 +48,7 @@ for which its arguments must be specified each time it is used. This is useful w
 PDEs for example, where one may need to use `u(t, x)` in the equations, but will
 need to be able to write `u(t, 0.0)` to define a boundary condition at `x = 0`.
 
-## Variable metadata [Experimental/TODO]
+## Variable metadata
 
 In many engineering systems, some variables act like “flows” while others do not.
 For example, in circuit models you have current which flows, and the related
@@ -69,15 +69,17 @@ the metadata. One can get and set metadata by
 ```julia
 julia> @variables x [unit = u"m^3/s"];
 
-julia> hasmetadata(x, Symbolics.option_to_metadata_type(Val(:unit)))
+julia> hasmetadata(x, VariableUnit)
 true
 
-julia> getmetadata(x, Symbolics.option_to_metadata_type(Val(:unit)))
+julia> ModelingToolkit.get_unit(x)
 m³ s⁻¹
 
-julia> x = setmetadata(x, Symbolics.option_to_metadata_type(Val(:unit)), u"m/s")
+julia> x = setmetadata(x, VariableUnit, u"m/s")
 x
 
-julia> getmetadata(x, Symbolics.option_to_metadata_type(Val(:unit)))
+julia> ModelingToolkit.get_unit(x)
 m s⁻¹
 ```
+
+See [Symbolic Metadata](@ref symbolic_metadata) for more details on variable metadata.
diff --git a/docs/src/basics/MTKModel_Connector.md b/docs/src/basics/MTKModel_Connector.md
index 8cc106b05d..9645260ed6 100644
--- a/docs/src/basics/MTKModel_Connector.md
+++ b/docs/src/basics/MTKModel_Connector.md
@@ -1,6 +1,23 @@
-# Defining components with `@mtkmodel`
+# [Components and Connectors](@id mtkmodel_connector)
 
-`@mtkmodel` is a convenience macro to define ModelingToolkit components. It returns `ModelingToolkit.Model`, which includes a constructor that returns an ODESystem, a `structure` dictionary with metadata and flag `isconnector` which is set to `false`.
+## MTK Model
+
+MTK represents components and connectors with `Model`.
+
+```@docs
+ModelingToolkit.Model
+```
+
+## Components
+
+Components are models from various domains. These models contain states and their
+equations.
+
+### [Defining components with `@mtkmodel`](@id mtkmodel)
+
+`@mtkmodel` is a convenience macro to define components. It returns
+`ModelingToolkit.Model`, which includes a constructor that returns the ODESystem, a
+`structure` dictionary with metadata, and flag `isconnector` which is set to `false`.
 
 ## What can an MTK-Model definition have?
 
@@ -56,9 +73,9 @@ end
 
   - Parameters and variables are declared with respective begin blocks.
   - Variables must be functions of an independent variable.
-  - Optionally, default values and metadata can be specified for these parameters and variables. See `ModelB` in the above example.
+  - Optionally, initial guess and metadata can be specified for these parameters and variables. See `ModelB` in the above example.
   - Along with creating parameters and variables, keyword arguments of same name with default value `nothing` are created.
-  - Whenever a parameter or variable has default value, for example `v(t) = 0.0`, a symbolic variable named `v` with default value 0.0 and a keyword argument `v`, with default value `nothing` are created. <br> This way, users can optionally pass new value of `v` while creating a component.
+  - Whenever a parameter or variable has initial value, for example `v(t) = 0.0`, a symbolic variable named `v` with initial value 0.0 and a keyword argument `v`, with default value `nothing` are created. <br> This way, users can optionally pass new value of `v` while creating a component.
 
 ```julia
 julia > @named model_c = ModelC(; v = 2.0);
@@ -72,7 +89,7 @@ julia > ModelingToolkit.getdefault(model_c.v)
   - This block is for non symbolic input arguements. These are for inputs that usually are not meant to be part of components; but influence how they are defined. One can list inputs like boolean flags, functions etc... here.
   - Whenever default values are specified, unlike parameters/variables, they are reflected in the keyword argument list.
 
-### `@extend` block
+#### `@extend` begin block
 
 To extend a partial system,
 
@@ -86,7 +103,8 @@ julia> @named model_c = ModelC(; p1 = 2.0)
 
 ```
 
-However, as `p2` isn't listed in the model definition, its default can't be modified by users.
+However, as `p2` isn't listed in the model definition, its initial guess can't
+specified while creating an instance of `ModelC`.
 
 ### `@components` begin block
 
@@ -110,13 +128,26 @@ And as `k2` isn't listed in the sub-component definition of `ModelC`, its defaul
 
   - Any other Julia operations can be included with dedicated begin blocks.
 
-## Defining connectors with `@connector`
+## Connectors
+
+Connectors are special models that can be used to connect different components together.
+MTK provides 3 distinct connectors:
+
+  - `DomainConnector`: A connector which has only one state which is of `Flow` type,
+    specified by `[connect = Flow]`.
+  - `StreamConnector`: A connector which has atleast one stream variable, specified by
+    `[connect = Stream]`. A `StreamConnector` must have exactly one flow variable.
+  - `RegularConnector`: Connectors that don't fall under above categories.
 
-`@connector` returns `ModelingToolkit.Model`. It includes a constructor that returns a connector ODESystem, a `structure` dictionary with metadata and flag `isconnector` which is set to `true`.
+### [Defining connectors with `@connector`](@id connector)
+
+`@connector` returns `ModelingToolkit.Model`. It includes a constructor that returns
+a connector ODESystem, a `structure` dictionary with metadata, and flag `isconnector`
+which is set to `true`.
 
 A simple connector can be defined with syntax similar to following example:
 
-```julia
+```@example connector
 using ModelingToolkit
 
 @connector Pin begin
@@ -125,16 +156,47 @@ using ModelingToolkit
 end
 ```
 
-  - Variables (as function of independent variable) are listed out in the definition. These variables can optionally have default values and metadata like `descrption`, `connect` and so on.
+Variables (as functions of independent variable) are listed out in the definition. These variables can optionally have initial values and metadata like `description`, `connect` and so on. For more details on setting metadata, check out [Symbolic Metadata](@ref symbolic_metadata).
 
-`@connector`s accepts begin blocks of `@components`, `@equations`, `@extend`, `@parameters`, `@structural_parameters`, `@variables`. These keywords mean the same as described above for `@mtkmodel`.
+Similar to `@mtkmodel`, `@connector` accepts begin blocks of `@components`, `@equations`, `@extend`, `@parameters`, `@structural_parameters`, `@variables`. These keywords mean the same as described above for `@mtkmodel`.
+For example, the following `HydraulicFluid` connector is defined with parameters, variables and equations.
+
+```@example connector
+@connector HydraulicFluid begin
+    @parameters begin
+        ρ = 997
+        β = 2.09e9
+        μ = 0.0010016
+        n = 1
+        let_gas = 1
+        ρ_gas = 0.0073955
+        p_gas = -1000
+    end
+    @variables begin
+        dm(t) = 0.0, [connect = Flow]
+    end
+    @equations begin
+        dm ~ 0
+    end
+end
+```
 
 !!! note
+    
     For more examples of usage, checkout [ModelingToolkitStandardLibrary.jl](https://github.com/SciML/ModelingToolkitStandardLibrary.jl/)
 
-### What's a `structure` dictionary?
+## More on `Model.structure`
+
+`structure` stores metadata that describes composition of a model. It includes:
 
-For components defined with `@mtkmodel` or `@connector`, a dictionary with metadata is created. It lists `:components` (sub-component list), `:extend` (the extended states and base system), `:parameters`, `:variables`, ``:kwargs`` (list of keyword arguments), `:independent_variable`, `:equations`.
+  - `:components`: List of sub-components in the form of [[name, sub_component_name],...].
+  - `:extend`: The list of extended states, name given to the base system, and name of the base system.
+  - `:structural_parameters`: Dictionary of structural parameters mapped to their default values.
+  - `:parameters`: Dictionary of symbolic parameters mapped to their metadata.
+  - `:variables`: Dictionary of symbolic variables mapped to their metadata.
+  - `:kwargs`: Dictionary of keyword arguments mapped to their default values.
+  - `:independent_variable`: Independent variable, which is added while generating the Model.
+  - `:equations`: List of equations (represented as strings).
 
 For example, the structure of `ModelC` is:
 
diff --git a/docs/src/basics/Variable_metadata.md b/docs/src/basics/Variable_metadata.md
index 53b8eb1fcf..d186cdf959 100644
--- a/docs/src/basics/Variable_metadata.md
+++ b/docs/src/basics/Variable_metadata.md
@@ -1,4 +1,4 @@
-# Symbolic metadata
+# [Symbolic Metadata](@id symbolic_metadata)
 
 It is possible to add metadata to symbolic variables, the metadata will be displayed when calling help on a variable.
 
@@ -39,6 +39,22 @@ help?> u
   Symbolics.VariableSource: (:variables, :u)
 ```
 
+## Connect
+
+Variables in connectors can have `connect` metadata which describes the type of connections.
+
+`Flow` is used for variables that represent physical quantities that "flow" ex:
+current in a resistor. These variables sum up to zero in connections.
+
+`Stream` can be specified for variables that flow bi-directionally.
+
+```@example connect
+using ModelingToolkit
+
+@variables t, i(t) [connect = Flow]
+@variables k(t) [connect = Stream]
+```
+
 ## Input or output
 
 Designate a variable as either an input or an output using the following
diff --git a/docs/src/index.md b/docs/src/index.md
index fc30b36714..1820bb40d6 100644
--- a/docs/src/index.md
+++ b/docs/src/index.md
@@ -173,13 +173,16 @@ system:
 
   - Please refer to the
     [SciML ColPrac: Contributor's Guide on Collaborative Practices for Community Packages](https://github.com/SciML/ColPrac/blob/master/README.md)
-    for guidance on PRs, issues, and other matters relating to contributing to ModelingToolkit.
+    for guidance on PRs, issues, and other matters relating to contributing to SciML.
 
+  - See the [SciML Style Guide](https://github.com/SciML/SciMLStyle) for common coding practices and other style decisions.
   - There are a few community forums:
     
-      + The #diffeq-bridged channel in the [Julia Slack](https://julialang.org/slack/)
-      + [JuliaDiffEq](https://gitter.im/JuliaDiffEq/Lobby) on Gitter
-      + On the Julia Discourse forums (look for the [modelingtoolkit tag](https://discourse.julialang.org/tag/modelingtoolkit))
+      + The #diffeq-bridged and #sciml-bridged channels in the
+        [Julia Slack](https://julialang.org/slack/)
+      + The #diffeq-bridged and #sciml-bridged channels in the
+        [Julia Zulip](https://julialang.zulipchat.com/#narrow/stream/279055-sciml-bridged)
+      + On the [Julia Discourse forums](https://discourse.julialang.org)
       + See also [SciML Community page](https://sciml.ai/community/)
 
 ## Reproducibility
diff --git a/docs/src/systems/DiscreteSystem.md b/docs/src/systems/DiscreteSystem.md
index c9d4b81c48..d9934d5fb7 100644
--- a/docs/src/systems/DiscreteSystem.md
+++ b/docs/src/systems/DiscreteSystem.md
@@ -23,7 +23,7 @@ DiscreteSystem
 
 ```@docs
 DiscreteFunction(sys::ModelingToolkit.DiscreteSystem, args...)
-DiscreteProblem(sys::ModelingToolkit.DiscreteSystem, args...)
+DiscreteProblem(sys::ModelingToolkit.DiscreteSystem, u0map, tspan)
 ```
 
 ## Expression Constructors
diff --git a/docs/src/systems/JumpSystem.md b/docs/src/systems/JumpSystem.md
index 3dbeb8500d..68f0de2a5c 100644
--- a/docs/src/systems/JumpSystem.md
+++ b/docs/src/systems/JumpSystem.md
@@ -15,7 +15,7 @@ JumpSystem
 
 ## Transformations
 
-```@docs
+```@docs; canonical=false
 structural_simplify
 ```
 
@@ -23,7 +23,9 @@ structural_simplify
 
 ## Problem Constructors
 
+```@docs; canonical=false
+DiscreteProblem(sys::ModelingToolkit.DiscreteSystem, u0map, tspan)
+```
 ```@docs
-SciMLBase.DiscreteProblem(sys::JumpSystem,args...)
-JumpProcesses.JumpProblem(sys::JumpSystem,args...)
+JumpProblem(sys::JumpSystem, prob, aggregator)
 ```
diff --git a/docs/src/systems/NonlinearSystem.md b/docs/src/systems/NonlinearSystem.md
index d14e6a035c..dfdb8db508 100644
--- a/docs/src/systems/NonlinearSystem.md
+++ b/docs/src/systems/NonlinearSystem.md
@@ -15,7 +15,7 @@ NonlinearSystem
 
 ## Transformations
 
-```@docs
+```@docs; canonical=false
 structural_simplify
 alias_elimination
 tearing
@@ -23,14 +23,14 @@ tearing
 
 ## Analyses
 
-```@docs
+```@docs; canonical=false
 ModelingToolkit.isaffine
 ModelingToolkit.islinear
 ```
 
 ## Applicable Calculation and Generation Functions
 
-```julia
+```@docs; canonical=false
 calculate_jacobian
 generate_jacobian
 jacobian_sparsity
diff --git a/docs/src/systems/ODESystem.md b/docs/src/systems/ODESystem.md
index 207ab098e7..f79263436f 100644
--- a/docs/src/systems/ODESystem.md
+++ b/docs/src/systems/ODESystem.md
@@ -35,7 +35,7 @@ ModelingToolkit.isaffine
 
 ## Applicable Calculation and Generation Functions
 
-```julia
+```@docs; canonical=false
 calculate_jacobian
 calculate_tgrad
 calculate_factorized_W
@@ -56,7 +56,7 @@ SteadyStateProblem(sys::ModelingToolkit.AbstractODESystem, args...)
 ## Torn Problem Constructors
 
 ```@docs
-ODAEProblem(sys::ModelingToolkit.AbstractODESystem, args...)
+ODAEProblem
 ```
 
 ## Expression Constructors
diff --git a/docs/src/systems/SDESystem.md b/docs/src/systems/SDESystem.md
index 57b89b304a..17f062c5ae 100644
--- a/docs/src/systems/SDESystem.md
+++ b/docs/src/systems/SDESystem.md
@@ -22,17 +22,19 @@ sde = SDESystem(ode, noiseeqs)
 
 ## Transformations
 
-```@docs
+```@docs; canonical=false
 structural_simplify
 alias_elimination
-Girsanov_transform
+```
+```@docs
+ModelingToolkit.Girsanov_transform
 ```
 
 ## Analyses
 
 ## Applicable Calculation and Generation Functions
 
-```julia
+```@docs; canonical=false
 calculate_jacobian
 calculate_tgrad
 calculate_factorized_W
diff --git a/docs/src/tutorials/acausal_components.md b/docs/src/tutorials/acausal_components.md
index cb2031d91a..508912c860 100644
--- a/docs/src/tutorials/acausal_components.md
+++ b/docs/src/tutorials/acausal_components.md
@@ -23,80 +23,89 @@ equalities before solving. Let's see this in action.
 using ModelingToolkit, Plots, DifferentialEquations
 
 @variables t
-@connector function Pin(; name)
-    sts = @variables v(t)=1.0 i(t)=1.0 [connect = Flow]
-    ODESystem(Equation[], t, sts, []; name = name)
+@connector Pin begin
+    v(t)
+    i(t), [connect = Flow]
 end
 
-@component function Ground(; name)
-    @named g = Pin()
-    eqs = [g.v ~ 0]
-    compose(ODESystem(eqs, t, [], []; name = name), g)
+@mtkmodel Ground begin
+    @components begin
+        g = Pin()
+    end
+    @equations begin
+        g.v ~ 0
+    end
 end
 
-@component function OnePort(; name)
-    @named p = Pin()
-    @named n = Pin()
-    sts = @variables v(t)=1.0 i(t)=1.0
-    eqs = [v ~ p.v - n.v
+@mtkmodel OnePort begin
+    @components begin
+        p = Pin()
+        n = Pin()
+    end
+    @variables begin
+        v(t)
+        i(t)
+    end
+    @equations begin
+        v ~ p.v - n.v
         0 ~ p.i + n.i
-        i ~ p.i]
-    compose(ODESystem(eqs, t, sts, []; name = name), p, n)
+        i ~ p.i
+    end
 end
 
-@component function Resistor(; name, R = 1.0)
-    @named oneport = OnePort()
-    @unpack v, i = oneport
-    ps = @parameters R = R
-    eqs = [
-        v ~ i * R,
-    ]
-    extend(ODESystem(eqs, t, [], ps; name = name), oneport)
+@mtkmodel Resistor begin
+    @extend OnePort()
+    @parameters begin
+        R = 1.0 # Sets the default resistance
+    end
+    @equations begin
+        v ~ i * R
+    end
 end
 
-@component function Capacitor(; name, C = 1.0)
-    @named oneport = OnePort()
-    @unpack v, i = oneport
-    ps = @parameters C = C
-    D = Differential(t)
-    eqs = [
-        D(v) ~ i / C,
-    ]
-    extend(ODESystem(eqs, t, [], ps; name = name), oneport)
+D = Differential(t)
+
+@mtkmodel Capacitor begin
+    @extend OnePort()
+    @parameters begin
+        C = 1.0
+    end
+    @equations begin
+        D(v) ~ i / C
+    end
+end
+
+@mtkmodel ConstantVoltage begin
+    @extend OnePort()
+    @parameters begin
+        V = 1.0
+    end
+    @equations begin
+        V ~ v
+    end
 end
 
-@component function ConstantVoltage(; name, V = 1.0)
-    @named oneport = OnePort()
-    @unpack v = oneport
-    ps = @parameters V = V
-    eqs = [
-        V ~ v,
-    ]
-    extend(ODESystem(eqs, t, [], ps; name = name), oneport)
+@mtkmodel RCModel begin
+    @components begin
+        resistor = Resistor(R = 1.0)
+        capacitor = Capacitor(C = 1.0)
+        source = ConstantVoltage(V = 1.0)
+        ground = Ground()
+    end
+    @equations begin
+        connect(source.p, resistor.p)
+        connect(resistor.n, capacitor.p)
+        connect(capacitor.n, source.n)
+        connect(capacitor.n, ground.g)
+    end
 end
 
-R = 1.0
-C = 1.0
-V = 1.0
-@named resistor = Resistor(R = R)
-@named capacitor = Capacitor(C = C)
-@named source = ConstantVoltage(V = V)
-@named ground = Ground()
-
-rc_eqs = [connect(source.p, resistor.p)
-    connect(resistor.n, capacitor.p)
-    connect(capacitor.n, source.n)
-    connect(capacitor.n, ground.g)]
-
-@named _rc_model = ODESystem(rc_eqs, t)
-@named rc_model = compose(_rc_model,
-    [resistor, capacitor, source, ground])
-sys = structural_simplify(rc_model)
+@mtkbuild rc_model = RCModel(resistor.R = 2.0)
 u0 = [
-    capacitor.v => 0.0,
+    rc_model.capacitor.v => 0.0,
 ]
-prob = ODAEProblem(sys, u0, (0, 10.0))
-sol = solve(prob, Tsit5())
+prob = ODEProblem(rc_model, u0, (0, 10.0))
+sol = solve(prob)
 plot(sol)
 ```
 
@@ -108,7 +117,7 @@ We wish to build the following RC circuit by building individual components and
 
 ### Building the Component Library
 
-For each of our components, we use a Julia function which emits an `ODESystem`.
+For each of our components, we use ModelingToolkit `Model` that emits an `ODESystem`.
 At the top, we start with defining the fundamental qualities of an electric
 circuit component. At every input and output pin, a circuit component has
 two values: the current at the pin and the voltage. Thus we define the `Pin`
@@ -119,40 +128,35 @@ i.e. that currents sum to zero and voltages across the pins are equal.
 default, variables are equal in a connection.
 
 ```@example acausal
-@connector function Pin(; name)
-    sts = @variables v(t)=1.0 i(t)=1.0 [connect = Flow]
-    ODESystem(Equation[], t, sts, []; name = name)
+@connector Pin begin
+    v(t)
+    i(t), [connect = Flow]
 end
 ```
 
 Note that this is an incompletely specified ODESystem: it cannot be simulated
 on its own because the equations for `v(t)` and `i(t)` are unknown. Instead,
 this just gives a common syntax for receiving this pair with some default
-values. Notice that in a component, we define the `name` as a keyword argument:
-this is because later we will generate different `Pin` objects with different
-names to correspond to duplicates of this topology with unique variables.
-One can then construct a `Pin` like:
-
-```@example acausal
-Pin(name = :mypin1)
-```
-
-or equivalently, using the `@named` helper macro:
+values.
+One can then construct a `Pin` using the `@named` helper macro:
 
 ```@example acausal
 @named mypin1 = Pin()
 ```
 
 Next, we build our ground node. A ground node is just a pin that is connected
-to a constant voltage reservoir, typically taken to be `V=0`. Thus to define
+to a constant voltage reservoir, typically taken to be `V = 0`. Thus to define
 this component, we generate an `ODESystem` with a `Pin` subcomponent and specify
 that the voltage in such a `Pin` is equal to zero. This gives:
 
 ```@example acausal
-@component function Ground(; name)
-    @named g = Pin()
-    eqs = [g.v ~ 0]
-    compose(ODESystem(eqs, t, [], []; name = name), g)
+@mtkmodel Ground begin
+    @components begin
+        g = Pin()
+    end
+    @equations begin
+        g.v ~ 0
+    end
 end
 ```
 
@@ -163,14 +167,20 @@ zero, and the current of the component equals to the current of the positive
 pin.
 
 ```@example acausal
-@component function OnePort(; name)
-    @named p = Pin()
-    @named n = Pin()
-    sts = @variables v(t)=1.0 i(t)=1.0
-    eqs = [v ~ p.v - n.v
+@mtkmodel OnePort begin
+    @components begin
+        p = Pin()
+        n = Pin()
+    end
+    @variables begin
+        v(t)
+        i(t)
+    end
+    @equations begin
+        v ~ p.v - n.v
         0 ~ p.i + n.i
-        i ~ p.i]
-    compose(ODESystem(eqs, t, sts, []; name = name), p, n)
+        i ~ p.i
+    end
 end
 ```
 
@@ -182,38 +192,37 @@ of charge we know that the current in must equal the current out, which means
 zero. This leads to our resistor equations:
 
 ```@example acausal
-@component function Resistor(; name, R = 1.0)
-    @named oneport = OnePort()
-    @unpack v, i = oneport
-    ps = @parameters R = R
-    eqs = [
-        v ~ i * R,
-    ]
-    extend(ODESystem(eqs, t, [], ps; name = name), oneport)
+@mtkmodel Resistor begin
+    @extend OnePort()
+    @parameters begin
+        R = 1.0
+    end
+    @equations begin
+        v ~ i * R
+    end
 end
 ```
 
 Notice that we have created this system with a default parameter `R` for the
 resistor's resistance. By doing so, if the resistance of this resistor is not
 overridden by a higher level default or overridden at `ODEProblem` construction
-time, this will be the value of the resistance. Also, note the use of `@unpack`
-and `extend`. For the `Resistor`, we want to simply inherit `OnePort`'s
-equations and states and extend them with a new equation. ModelingToolkit makes
-a new namespaced variable `oneport₊v(t)` when using the syntax `oneport.v`, and
-we can use `@unpack` to avoid the namespacing.
+time, this will be the value of the resistance. Also, note the use of `@extend`.
+For the `Resistor`, we want to simply inherit `OnePort`'s
+equations and states and extend them with a new equation. Note that `v`, `i` are not namespaced as `oneport.v` or `oneport.i`.
 
 Using our knowledge of circuits, we similarly construct the `Capacitor`:
 
 ```@example acausal
-@component function Capacitor(; name, C = 1.0)
-    @named oneport = OnePort()
-    @unpack v, i = oneport
-    ps = @parameters C = C
-    D = Differential(t)
-    eqs = [
-        D(v) ~ i / C,
-    ]
-    extend(ODESystem(eqs, t, [], ps; name = name), oneport)
+D = Differential(t)
+
+@mtkmodel Capacitor begin
+    @extend OnePort()
+    @parameters begin
+        C = 1.0
+    end
+    @equations begin
+        D(v) ~ i / C
+    end
 end
 ```
 
@@ -223,55 +232,52 @@ constant voltage, essentially generating the electric current. We would then
 model this as:
 
 ```@example acausal
-@component function ConstantVoltage(; name, V = 1.0)
-    @named oneport = OnePort()
-    @unpack v = oneport
-    ps = @parameters V = V
-    eqs = [
-        V ~ v,
-    ]
-    extend(ODESystem(eqs, t, [], ps; name = name), oneport)
+@mtkmodel ConstantVoltage begin
+    @extend OnePort()
+    @parameters begin
+        V = 1.0
+    end
+    @equations begin
+        V ~ v
+    end
 end
 ```
 
+Note that as we are extending only `v` from `OnePort`, it is explicitly specified as a tuple.
+
 ### Connecting and Simulating Our Electric Circuit
 
 Now we are ready to simulate our circuit. Let's build our four components:
 a `resistor`, `capacitor`, `source`, and `ground` term. For simplicity, we will
-make all of our parameter values 1. This is done by:
-
-```@example acausal
-R = 1.0
-C = 1.0
-V = 1.0
-@named resistor = Resistor(R = R)
-@named capacitor = Capacitor(C = C)
-@named source = ConstantVoltage(V = V)
-@named ground = Ground()
-```
-
-Subsequently, we will connect the pieces of our circuit together. Let's connect the
-positive pin of the resistor to the source, the negative pin of the resistor
-to the capacitor, and the negative pin of the capacitor to a junction between
-the source and the ground. This would mean our connection equations are:
+make all of our parameter values 1.0. As `resistor`, `capacitor`, `source` lists
+`R`, `C`, `V` in their argument list, they are promoted as arguments of RCModel as
+`resistor.R`, `capacitor.C`, `source.V`
 
 ```@example acausal
-rc_eqs = [connect(source.p, resistor.p)
-    connect(resistor.n, capacitor.p)
-    connect(capacitor.n, source.n)
-    connect(capacitor.n, ground.g)]
+@mtkmodel RCModel begin
+    @components begin
+        resistor = Resistor(R = 1.0)
+        capacitor = Capacitor(C = 1.0)
+        source = ConstantVoltage(V = 1.0)
+        ground = Ground()
+    end
+    @equations begin
+        connect(source.p, resistor.p)
+        connect(resistor.n, capacitor.p)
+        connect(capacitor.n, source.n)
+        connect(capacitor.n, ground.g)
+    end
+end
 ```
 
-Finally, we build our four-component model with these connection rules:
+We can create a RCModel component with `@named`. And using `subcomponent_name.parameter` we can set
+the parameters or defaults values of variables of subcomponents.
 
 ```@example acausal
-@named _rc_model = ODESystem(rc_eqs, t)
-@named rc_model = compose(_rc_model,
-    [resistor, capacitor, source, ground])
+@mtkbuild rc_model = RCModel(resistor.R = 2.0)
 ```
 
-Note that we can also specify the subsystems in a vector. This model is acausal
-because we have not specified anything about the causality of the model. We have
+This model is acausal because we have not specified anything about the causality of the model. We have
 simply specified what is true about each of the variables. This forms a system
 of differential-algebraic equations (DAEs) which define the evolution of each
 state of the system. The equations are:
@@ -292,26 +298,15 @@ and the parameters are:
 parameters(rc_model)
 ```
 
-## Simplifying and Solving this System
-
-This system could be solved directly as a DAE using [one of the DAE solvers
-from DifferentialEquations.jl](https://docs.sciml.ai/DiffEqDocs/stable/solvers/dae_solve/).
-However, let's take a second to symbolically simplify the system before doing the
-solve. Although we can use ODE solvers that handle mass matrices to solve the
-above system directly, we want to run the `structural_simplify` function first,
-as it eliminates many unnecessary variables to build the leanest numerical
-representation of the system. Let's see what it does here:
+The observed equations are:
 
 ```@example acausal
-sys = structural_simplify(rc_model)
-equations(sys)
+observed(rc_model)
 ```
 
-```@example acausal
-states(sys)
-```
+## Solving this System
 
-After structural simplification, we are left with a system of only two equations
+We are left with a system of only two equations
 with two state variables. One of the equations is a differential equation,
 while the other is an algebraic equation. We can then give the values for the
 initial conditions of our states, and solve the system by converting it to
@@ -320,35 +315,21 @@ DAE solver](https://docs.sciml.ai/DiffEqDocs/stable/solvers/dae_solve/#OrdinaryD
 This is done as follows:
 
 ```@example acausal
-u0 = [capacitor.v => 0.0
-    capacitor.p.i => 0.0]
-prob = ODEProblem(sys, u0, (0, 10.0))
-sol = solve(prob, Rodas4())
-plot(sol)
-```
+u0 = [rc_model.capacitor.v => 0.0
+    rc_model.capacitor.p.i => 0.0]
 
-Since we have run `structural_simplify`, MTK can numerically solve all the
-unreduced algebraic equations using the `ODAEProblem` (note the
-letter `A`):
-
-```@example acausal
-u0 = [
-    capacitor.v => 0.0,
-]
-prob = ODAEProblem(sys, u0, (0, 10.0))
-sol = solve(prob, Rodas4())
+prob = ODEProblem(rc_model, u0, (0, 10.0))
+sol = solve(prob)
 plot(sol)
 ```
 
-Notice that this solves the whole system by only solving for one variable!
-
 However, what if we wanted to plot the timeseries of a different variable? Do
 not worry, that information was not thrown away! Instead, transformations
-like `structural_simplify` simply change state variables into `observed`
-variables. Let's see what our observed variables are:
+like `structural_simplify` simply change state variables into observables which are
+defined by `observed` equations.
 
 ```@example acausal
-observed(sys)
+observed(rc_model)
 ```
 
 These are explicit algebraic equations which can then be used to reconstruct
@@ -360,11 +341,11 @@ The solution object can be accessed via its symbols. For example, let's retrieve
 the voltage of the resistor over time:
 
 ```@example acausal
-sol[resistor.v]
+sol[rc_model.resistor.v]
 ```
 
 or we can plot the timeseries of the resistor's voltage:
 
 ```@example acausal
-plot(sol, idxs = [resistor.v])
+plot(sol, idxs = [rc_model.resistor.v])
 ```
diff --git a/docs/src/tutorials/bifurcation_diagram_computation.md b/docs/src/tutorials/bifurcation_diagram_computation.md
new file mode 100644
index 0000000000..1bb6e33935
--- /dev/null
+++ b/docs/src/tutorials/bifurcation_diagram_computation.md
@@ -0,0 +1,93 @@
+# [Bifurcation Diagrams](@id bifurcation_diagrams)
+Bifurcation diagrams describes how, for a dynamic system, the quantity and quality of its steady states changes with a parameter's value. These can be computed through the [BifurcationKit.jl](https://github.com/bifurcationkit/BifurcationKit.jl) package. ModelingToolkit provides a simple interface for creating BifurcationKit compatible `BifurcationProblem`s from `NonlinearSystem`s and `ODESystem`s. All teh features provided by BifurcationKit can then be applied to these systems. This tutorial provides a brief introduction for these features, with BifurcationKit.jl providing [a more extensive documentation](https://bifurcationkit.github.io/BifurcationKitDocs.jl/stable/).
+
+### Creating a `BifurcationProblem`
+Let us first consider a simple `NonlinearSystem`:
+```@example Bif1
+using ModelingToolkit
+@variables t x(t) y(t)
+@parameters μ α
+eqs = [0 ~ μ*x - x^3 + α*y,
+        0 ~ -y]
+@named nsys = NonlinearSystem(eqs, [x, y], [μ, α])
+```
+we wish to compute a bifurcation diagram for this system as we vary the parameter `μ`. For this, we need to provide the following information:
+1. The system for which we wish to compute the bifurcation diagram (`nsys`).
+2. The parameter which we wish to vary (`μ`).
+3. The parameter set for which we want to compute the bifurcation diagram.
+4. An initial guess of the state of the system for which there is a steady state at our provided parameter value.
+5. The variable which value we wish to plot in the bifurcation diagram (this argument is optional, if not provided, BifurcationKit default plot functions are used).
+
+We declare this additional information:
+```@example Bif1
+bif_par = μ
+p_start = [μ => -1.0, α => 1.0]
+u0_guess = [x => 1.0, y => 1.0]
+plot_var = x;
+```
+For the initial state guess (`u0_guess`), typically any value can be provided, however, read BifurcatioKit's documentation for more details.
+
+We can now create our `BifurcationProblem`, which can be provided as input to BifurcationKit's various functions.
+```@example Bif1
+using BifurcationKit
+bprob = BifurcationProblem(nsys, u0_guess, p_start, bif_par; plot_var=plot_var, jac=false)
+```
+Here, the `jac` argument (by default set to `true`) sets whenever to provide BifurcationKit with a Jacobian or not.
+
+
+### Computing a bifurcation diagram
+
+Let us consider the `BifurcationProblem` from the last section. If we wish to compute the corresponding bifurcation diagram we must first declare various settings used by BifurcationKit to compute the diagram. These are stored in a `ContinuationPar` structure (which also contain a `NewtonPar` structure). 
+```@example Bif1
+p_span = (-4.0, 6.0)
+opt_newton = NewtonPar(tol = 1e-9, max_iterations = 20)
+opts_br = ContinuationPar(dsmin = 0.001, dsmax = 0.05, ds = 0.01,
+	max_steps = 100, nev = 2, newton_options = opt_newton,
+	p_min = p_span[1], p_max = p_span[2],
+	detect_bifurcation = 3, n_inversion = 4, tol_bisection_eigenvalue = 1e-8, dsmin_bisection = 1e-9);
+```
+Here, `p_span` sets the interval over which we wish to compute the diagram.
+
+Next, we can use this as input to our bifurcation diagram, and then plot it.
+```@example Bif1
+bf = bifurcationdiagram(bprob, PALC(), 2, (args...) -> opts_br; bothside=true)
+```
+Here, the value `2` sets how sub-branches of the diagram that BifurcationKit should compute. Generally, for bifurcation diagrams, it is recommended to use the `bothside=true` argument.
+```@example Bif1
+using Plots
+plot(bf; putspecialptlegend=false, markersize=2, plotfold=false, xguide="μ", yguide = "x")
+```
+Here, the system exhibits a pitchfork bifurcation at *μ=0.0*.
+
+### Using `ODESystem` inputs
+It is also possible to use `ODESystem`s (rather than `NonlinearSystem`s) as input to `BifurcationProblem`. Here follows a brief such example.
+
+```@example Bif2
+using BifurcationKit, ModelingToolkit, Plots
+
+@variables t x(t) y(t)
+@parameters μ
+D = Differential(t)
+eqs = [D(x) ~ μ*x - y - x*(x^2+y^2), 
+        D(y) ~ x + μ*y - y*(x^2+y^2)]
+@named osys = ODESystem(eqs, t)
+
+bif_par = μ
+plot_var = x
+p_start = [μ => 1.0]
+u0_guess = [x => 0.0, y=> 0.0]
+
+bprob = BifurcationProblem(osys, u0_guess, p_start, bif_par; plot_var=plot_var, jac=false)
+
+p_span = (-3.0, 3.0)
+opt_newton = NewtonPar(tol = 1e-9, max_iterations = 20)
+opts_br = ContinuationPar(dsmin = 0.001, dsmax = 0.05, ds = 0.01,
+	max_steps = 100, nev = 2, newton_options = opt_newton,
+	p_max = p_span[2], p_min = p_span[1],
+	detect_bifurcation = 3, n_inversion = 4, tol_bisection_eigenvalue = 1e-8, dsmin_bisection = 1e-9)
+
+bf = bifurcationdiagram(bprob, PALC(), 2, (args...) -> opts_br; bothside=true)
+using Plots
+plot(bf; putspecialptlegend=false, markersize=2, plotfold=false, xguide="μ", yguide = "x")
+```
+Here, the value of `x` in the steady state does not change, however, at `μ=0` a Hopf bifurcation occur and the steady state turn unstable.
\ No newline at end of file
diff --git a/docs/src/tutorials/ode_modeling.md b/docs/src/tutorials/ode_modeling.md
index 5763eeabf6..ae56d4f44d 100644
--- a/docs/src/tutorials/ode_modeling.md
+++ b/docs/src/tutorials/ode_modeling.md
@@ -1,29 +1,44 @@
 # Getting Started with ModelingToolkit.jl
 
-This is an introductory tutorial for ModelingToolkit (MTK).
-Some examples of Ordinary Differential Equations (ODE) are used to
-illustrate the basic user-facing functionality.
+This is an introductory tutorial for ModelingToolkit (MTK). We will demonstrate
+the basics of the package by demonstrating how to define and simulate simple
+Ordinary Differential Equation (ODE) systems.
+
+## Installing ModelingToolkit
+
+To install ModelingToolkit, use the Julia package manager. This can be done as follows:
+
+```julia
+using Pkg
+Pkg.add("ModelingToolkit")
+```
 
 ## Copy-Pastable Simplified Example
 
 A much deeper tutorial with forcing functions and sparse Jacobians is below.
 But if you want to just see some code and run, here's an example:
 
-```@example ode
+```@example first-mtkmodel
 using ModelingToolkit
 
-@variables t x(t)   # independent and dependent variables
-@parameters τ       # parameters
-@constants h = 1    # constants have an assigned value
-D = Differential(t) # define an operator for the differentiation w.r.t. time
-
-# your first ODE, consisting of a single equation, the equality indicated by ~
-@named fol = ODESystem([D(x) ~ (h - x) / τ])
+@variables t
+D = Differential(t)
+
+@mtkmodel FOL begin
+    @parameters begin
+        τ # parameters
+    end
+    @variables begin
+        x(t) # dependent variables
+    end
+    @equations begin
+        D(x) ~ (1 - x) / τ
+    end
+end
 
 using DifferentialEquations: solve
-
-prob = ODEProblem(fol, [x => 0.0], (0.0, 10.0), [τ => 3.0])
-# parameter `τ` can be assigned a value, but constant `h` cannot
+@mtkbuild fol = FOL()
+prob = ODEProblem(fol, [fol.x => 0.0], (0.0, 10.0), [fol.τ => 3.0])
 sol = solve(prob)
 
 using Plots
@@ -43,24 +58,35 @@ first-order lag element:
 
 Here, ``t`` is the independent variable (time), ``x(t)`` is the (scalar) state
 variable, ``f(t)`` is an external forcing function, and ``\tau`` is a
-parameter. In MTK, this system can be modelled as follows. For simplicity, we
-first set the forcing function to a time-independent value.
+parameter.
+In MTK, this system can be modelled as follows. For simplicity, we
+first set the forcing function to a time-independent value ``1``. And the
+independent variable ``t`` is automatically added by ``@mtkmodel``.
 
 ```@example ode2
 using ModelingToolkit
 
-@variables t x(t)  # independent and dependent variables
-@parameters τ       # parameters
-@constants h = 1    # constants
-D = Differential(t) # define an operator for the differentiation w.r.t. time
+@variables t
+D = Differential(t)
+
+@mtkmodel FOL begin
+    @parameters begin
+        τ # parameters
+    end
+    @variables begin
+        x(t) # dependent variables
+    end
+    @equations begin
+        D(x) ~ (1 - x) / τ
+    end
+end
 
-# your first ODE, consisting of a single equation, indicated by ~
-@named fol_model = ODESystem(D(x) ~ (h - x) / τ)
+@mtkbuild fol = FOL()
 ```
 
 Note that equations in MTK use the tilde character (`~`) as equality sign.
-Also note that the `@named` macro simply ensures that the symbolic name
-matches the name in the REPL. If omitted, you can directly set the `name` keyword.
+
+`@mtkbuild` creates an instance of `FOL` named as `fol`.
 
 After construction of the ODE, you can solve it using [DifferentialEquations.jl](https://docs.sciml.ai/DiffEqDocs/stable/):
 
@@ -68,7 +94,7 @@ After construction of the ODE, you can solve it using [DifferentialEquations.jl]
 using DifferentialEquations
 using Plots
 
-prob = ODEProblem(fol_model, [x => 0.0], (0.0, 10.0), [τ => 3.0])
+prob = ODEProblem(fol, [fol.x => 0.0], (0.0, 10.0), [fol.τ => 3.0])
 plot(solve(prob))
 ```
 
@@ -82,53 +108,80 @@ You could separate the calculation of the right-hand side, by introducing an
 intermediate variable `RHS`:
 
 ```@example ode2
-@variables RHS(t)
-@named fol_separate = ODESystem([RHS ~ (h - x) / τ,
-    D(x) ~ RHS])
+using ModelingToolkit
+
+@mtkmodel FOL begin
+    @parameters begin
+        τ # parameters
+    end
+    @variables begin
+        x(t) # dependent variables
+        RHS(t)
+    end
+    begin
+        D = Differential(t)
+    end
+    @equations begin
+        RHS ~ (1 - x) / τ
+        D(x) ~ RHS
+    end
+end
+
+@mtkbuild fol = FOL()
 ```
 
-To directly solve this system, you would have to create a Differential-Algebraic
-Equation (DAE) problem, since besides the differential equation, there is an
-additional algebraic equation now. However, this DAE system can obviously be
-transformed into the single ODE we used in the first example above. MTK achieves
-this by structural simplification:
+You can look at the equations by using the command `equations`:
 
 ```@example ode2
-fol_simplified = structural_simplify(fol_separate)
-equations(fol_simplified)
+equations(fol)
 ```
 
+Notice that there is only one equation in this system, `Differential(t)(x(t)) ~ RHS(t)`.
+The other equation was removed from the system and was transformed into an `observed`
+variable. Observed equations are variables which can be computed on-demand but are not
+necessary for the solution of the system, and thus MTK tracks it separately. One can
+check the observed equations via the `observed` function:
+
 ```@example ode2
-equations(fol_simplified) == equations(fol_model)
+observed(fol)
 ```
 
-You can extract the equations from a system using `equations` (and, in the same
-way, `states` and `parameters`). The simplified equation is exactly the same
-as the original one, so the simulation performance will also be the same.
-However, there is one difference. MTK does keep track of the eliminated
-algebraic variables as "observables" (see
-[Observables and Variable Elimination](@ref)).
-That means, MTK still knows how to calculate them out of the information available
+For more information on this process, see [Observables and Variable Elimination](@ref).
+
+MTK still knows how to calculate them out of the information available
 in a simulation result. The intermediate variable `RHS` therefore can be plotted
-along with the state variable. Note that this has to be requested explicitly,
-through:
+along with the state variable. Note that this has to be requested explicitly
+like as follows:
 
 ```@example ode2
-prob = ODEProblem(fol_simplified, [x => 0.0], (0.0, 10.0), [τ => 3.0])
+prob = ODEProblem(fol,
+    [fol.x => 0.0],
+    (0.0, 10.0),
+    [fol.τ => 3.0])
 sol = solve(prob)
-plot(sol, vars = [x, RHS])
+plot(sol, vars = [fol.x, fol.RHS])
 ```
 
-By default, `structural_simplify` also replaces symbolic `constants` with
-their default values. This allows additional simplifications not possible
-when using `parameters` (e.g., solution of linear equations by dividing out
-the constant's value, which cannot be done for parameters, since they may
-be zero).
+## Named Indexing of Solutions
+
+Note that the indexing of the solution similarly works via the names, and so to get
+the time series for `x`, one would do:
 
-Note that the indexing of the solution similarly works via the names, and so
-`sol[x]` gives the time-series for `x`, `sol[x,2:10]` gives the 2nd through 10th
-values of `x` matching `sol.t`, etc. Note that this works even for variables
-which have been eliminated, and thus `sol[RHS]` retrieves the values of `RHS`.
+```@example ode2
+sol[fol.x]
+```
+
+or to get the second value in the time series for `x`:
+
+```@example ode2
+sol[fol.x, 2]
+```
+
+Similarly, the time series for `RHS` can be retrieved using the same indexing:
+
+```@example ode2
+sol[fol.RHS]
+```
 
 ## Specifying a time-variable forcing function
 
@@ -136,8 +189,24 @@ What if the forcing function (the “external input”) ``f(t)`` is not constant
 Obviously, one could use an explicit, symbolic function of time:
 
 ```@example ode2
-@variables f(t)
-@named fol_variable_f = ODESystem([f ~ sin(t), D(x) ~ (f - x) / τ])
+@mtkmodel FOL begin
+    @parameters begin
+        τ # parameters
+    end
+    @variables begin
+        x(t) # dependent variables
+        f(t)
+    end
+    begin
+        D = Differential(t)
+    end
+    @equations begin
+        f ~ sin(t)
+        D(x) ~ (f - x) / τ
+    end
+end
+
+@named fol_variable_f = FOL()
 ```
 
 But often this function might not be available in an explicit form.
@@ -154,11 +223,34 @@ value_vector = randn(10)
 f_fun(t) = t >= 10 ? value_vector[end] : value_vector[Int(floor(t)) + 1]
 @register_symbolic f_fun(t)
 
-@named fol_external_f = ODESystem([f ~ f_fun(t), D(x) ~ (f - x) / τ])
-prob = ODEProblem(structural_simplify(fol_external_f), [x => 0.0], (0.0, 10.0), [τ => 0.75])
+@mtkmodel FOLExternalFunction begin
+    @parameters begin
+        τ # parameters
+    end
+    @variables begin
+        x(t) # dependent variables
+        f(t)
+    end
+    @structural_parameters begin
+        h = 1
+    end
+    begin
+        D = Differential(t)
+    end
+    @equations begin
+        f ~ f_fun(t)
+        D(x) ~ (f - x) / τ
+    end
+end
+
+@mtkbuild fol_external_f = FOLExternalFunction()
+prob = ODEProblem(fol_external_f,
+    [fol_external_f.x => 0.0],
+    (0.0, 10.0),
+    [fol_external_f.τ => 0.75])
 
 sol = solve(prob)
-plot(sol, vars = [x, f])
+plot(sol, vars = [fol_external_f.x, fol_external_f.f])
 ```
 
 ## Building component-based, hierarchical models
@@ -235,19 +327,43 @@ plot(solve(prob))
 
 More on this topic may be found in [Composing Models and Building Reusable Components](@ref acausal).
 
-## Defaults
+## Inital Guess
 
 It is often a good idea to specify reasonable values for the initial state and the
 parameters of a model component. Then, these do not have to be explicitly specified when constructing the `ODEProblem`.
 
 ```@example ode2
-function unitstep_fol_factory(; name)
+@mtkmodel UnitstepFOLFactory begin
+    @parameters begin
+        τ = 1.0
+    end
+    @variables begin
+        x(t) = 0.0
+    end
+    @equations begin
+        D(x) ~ (1 - x) / τ
+    end
+end
+```
+
+While defining the model `UnitstepFOLFactory`, an initial guess of 0.0 is assigned to `x(t)` and 1.0 to `τ`.
+Additionaly, these initial guesses can be modified while creating instances of `UnitstepFOLFactory` by passing arguements.
+
+```@example ode2
+@named fol = UnitstepFOLFactory(; x = 0.1)
+sol = ODEProblem(fol, [], (0.0, 5.0), []) |> solve
+```
+
+In non-DSL definitions, one can pass `defaults` dictionary to set the initial guess of the symbolic variables.
+
+```@example ode3
+using ModelingToolkit
+
+function UnitstepFOLFactory(; name)
     @parameters τ
     @variables t x(t)
     ODESystem(D(x) ~ (1 - x) / τ; name, defaults = Dict(x => 0.0, τ => 1.0))
 end
-
-ODEProblem(unitstep_fol_factory(name = :fol), [], (0.0, 5.0), []) |> solve
 ```
 
 Note that the defaults can be functions of the other variables, which is then
@@ -300,14 +416,9 @@ using the structural information. For more information, see the
 
 Here are some notes that may be helpful during your initial steps with MTK:
 
-  - Sometimes, the symbolic engine within MTK cannot correctly identify the
-    independent variable (e.g. time) out of all variables. In such a case, you
-    usually get an error that some variable(s) is "missing from variable map". In
-    most cases, it is then sufficient to specify the independent variable as second
-    argument to `ODESystem`, e.g. `ODESystem(eqs, t)`.
-  - A completely macro-free usage of MTK is possible and is discussed in a
-    separate tutorial. This is for package developers, since the macros are only
-    essential for automatic symbolic naming for modelers.
+  - The `@mtkmodel` macro is for high-level usage of MTK. However, in many cases you
+    may need to programmatically generate `ODESystem`s. If that's the case, check out
+    the [Programmatically Generating and Scripting ODESystems Tutorial](@ref programmatically).
   - Vector-valued parameters and variables are possible. A cleaner, more
     consistent treatment of these is still a work in progress, however. Once finished,
     this introductory tutorial will also cover this feature.
@@ -316,6 +427,8 @@ Where to go next?
 
   - Not sure how MTK relates to similar tools and packages? Read
     [Comparison of ModelingToolkit vs Equation-Based and Block Modeling Languages](@ref).
+  - For a more detailed explanation of `@mtkmodel` checkout
+    [Defining components with `@mtkmodel` and connectors with `@connectors`](@ref mtkmodel_connector)
   - Depending on what you want to do with MTK, have a look at some of the other
     **Symbolic Modeling Tutorials**.
   - If you want to automatically convert an existing function to a symbolic
diff --git a/docs/src/tutorials/optimization.md b/docs/src/tutorials/optimization.md
index 24b74c737b..b4ee03141b 100644
--- a/docs/src/tutorials/optimization.md
+++ b/docs/src/tutorials/optimization.md
@@ -69,8 +69,8 @@ cons = [
     x^2 + y^2 ≲ 1,
 ]
 @named sys = OptimizationSystem(loss, [x, y], [a, b], constraints = cons)
-u0 = [x => 1.0
-    y => 2.0]
+u0 = [x => 0.14
+    y => 0.14]
 prob = OptimizationProblem(sys, u0, grad = true, hess = true, cons_j = true, cons_h = true)
 solve(prob, IPNewton())
 ```
diff --git a/docs/src/tutorials/parameter_identifiability.md b/docs/src/tutorials/parameter_identifiability.md
index 670c3397cc..eb0837360b 100644
--- a/docs/src/tutorials/parameter_identifiability.md
+++ b/docs/src/tutorials/parameter_identifiability.md
@@ -31,24 +31,40 @@ We first define the parameters, variables, differential equations and the output
 ```@example SI
 using StructuralIdentifiability, ModelingToolkit
 
-# define parameters and variables
-@variables t x4(t) x5(t) x6(t) x7(t) y1(t) y2(t)
-@parameters k5 k6 k7 k8 k9 k10
+@variables t
 D = Differential(t)
 
-# define equations
-eqs = [
-    D(x4) ~ -k5 * x4 / (k6 + x4),
-    D(x5) ~ k5 * x4 / (k6 + x4) - k7 * x5 / (k8 + x5 + x6),
-    D(x6) ~ k7 * x5 / (k8 + x5 + x6) - k9 * x6 * (k10 - x6) / k10,
-    D(x7) ~ k9 * x6 * (k10 - x6) / k10,
-]
-
-# define the output functions (quantities that can be measured)
-measured_quantities = [y1 ~ x4, y2 ~ x5]
+@mtkmodel Biohydrogenation begin
+    @variables begin
+        x4(t)
+        x5(t)
+        x6(t)
+        x7(t)
+        y1(t), [output = true]
+        y2(t), [output = true]
+    end
+    @parameters begin
+        k5
+        k6
+        k7
+        k8
+        k9
+        k10
+    end
+    # define equations
+    @equations begin
+        D(x4) ~ -k5 * x4 / (k6 + x4)
+        D(x5) ~ k5 * x4 / (k6 + x4) - k7 * x5 / (k8 + x5 + x6)
+        D(x6) ~ k7 * x5 / (k8 + x5 + x6) - k9 * x6 * (k10 - x6) / k10
+        D(x7) ~ k9 * x6 * (k10 - x6) / k10
+        y1 ~ x4
+        y2 ~ x5
+    end
+end
 
 # define the system
-de = ODESystem(eqs, t, name = :Biohydrogenation)
+@named de = Biohydrogenation()
+de = complete(de)
 ```
 
 After that, we are ready to check the system for local identifiability:
@@ -56,8 +72,7 @@ After that, we are ready to check the system for local identifiability:
 ```@example SI
 # query local identifiability
 # we pass the ode-system
-local_id_all = assess_local_identifiability(de, measured_quantities = measured_quantities,
-    p = 0.99)
+local_id_all = assess_local_identifiability(de, p = 0.99)
 ```
 
 We can see that all states (except $x_7$) and all parameters are locally identifiable with probability 0.99.
@@ -65,9 +80,8 @@ We can see that all states (except $x_7$) and all parameters are locally identif
 Let's try to check specific parameters and their combinations
 
 ```@example SI
-to_check = [k5, k7, k10 / k9, k5 + k6]
-local_id_some = assess_local_identifiability(de, measured_quantities = measured_quantities,
-    funcs_to_check = to_check, p = 0.99)
+to_check = [de.k5, de.k7, de.k10 / de.k9, de.k5 + de.k6]
+local_id_some = assess_local_identifiability(de, funcs_to_check = to_check, p = 0.99)
 ```
 
 Notice that in this case, everything (except the state variable $x_7$) is locally identifiable, including combinations such as $k_{10}/k_9, k_5+k_6$
@@ -92,26 +106,43 @@ We will run a global identifiability check on this enzyme dynamics[^3] model. We
 
 Global identifiability needs information about local identifiability first, but the function we chose here will take care of that extra step for us.
 
-__Note__: as of writing this tutorial, UTF-symbols such as Greek characters are not supported by one of the project's dependencies, see [this issue](https://github.com/SciML/StructuralIdentifiability.jl/issues/43).
-
 ```@example SI2
 using StructuralIdentifiability, ModelingToolkit
-@parameters b c a beta g delta sigma
-@variables t x1(t) x2(t) x3(t) x4(t) y(t) y2(t)
-D = Differential(t)
-
-eqs = [
-    D(x1) ~ -b * x1 + 1 / (c + x4),
-    D(x2) ~ a * x1 - beta * x2,
-    D(x3) ~ g * x2 - delta * x3,
-    D(x4) ~ sigma * x4 * (g * x2 - delta * x3) / x3,
-]
 
-measured_quantities = [y ~ x1 + x2, y2 ~ x2]
-
-ode = ODESystem(eqs, t, name = :GoodwinOsc)
+@variables t
+D = Differential(t)
 
-global_id = assess_identifiability(ode, measured_quantities = measured_quantities)
+@mtkmodel GoodwinOsc begin
+    @parameters begin
+        b
+        c
+        α
+        β
+        γ
+        δ
+        σ
+    end
+    @variables begin
+        x1(t)
+        x2(t)
+        x3(t)
+        x4(t)
+        y(t), [output = true]
+        y2(t), [output = true]
+    end
+    @equations begin
+        D(x1) ~ -b * x1 + 1 / (c + x4)
+        D(x2) ~ α * x1 - β * x2
+        D(x3) ~ γ * x2 - δ * x3
+        D(x4) ~ σ * x4 * (γ * x2 - δ * x3) / x3
+        y ~ x1 + x2
+        y2 ~ x2
+    end
+end
+
+@named ode = GoodwinOsc()
+
+global_id = assess_identifiability(ode)
 ```
 
 We can see that only parameters `a, g` are unidentifiable, and everything else can be uniquely recovered.
@@ -120,25 +151,47 @@ Let us consider the same system but with two inputs, and we will find out identi
 
 ```@example SI3
 using StructuralIdentifiability, ModelingToolkit
-@parameters b c a beta g delta sigma
-@variables t x1(t) x2(t) x3(t) x4(t) y(t) y2(t) u1(t) [input = true] u2(t) [input = true]
+
+@variables t
 D = Differential(t)
 
-eqs = [
-    D(x1) ~ -b * x1 + 1 / (c + x4),
-    D(x2) ~ a * x1 - beta * x2 - u1,
-    D(x3) ~ g * x2 - delta * x3 + u2,
-    D(x4) ~ sigma * x4 * (g * x2 - delta * x3) / x3,
-]
-measured_quantities = [y ~ x1 + x2, y2 ~ x2]
+@mtkmodel GoodwinOscillator begin
+    @parameters begin
+        b
+        c
+        α
+        β
+        γ
+        δ
+        σ
+    end
+    @variables begin
+        x1(t)
+        x2(t)
+        x3(t)
+        x4(t)
+        y(t), [output = true]
+        y2(t), [output = true]
+        u1(t), [input = true]
+        u2(t), [input = true]
+    end
+    @equations begin
+        D(x1) ~ -b * x1 + 1 / (c + x4)
+        D(x2) ~ α * x1 - β * x2 - u1
+        D(x3) ~ γ * x2 - δ * x3 + u2
+        D(x4) ~ σ * x4 * (γ * x2 - δ * x3) / x3
+        y ~ x1 + x2
+        y2 ~ x2
+    end
+end
+
+@named ode = GoodwinOscillator()
+ode = complete(ode)
 
 # check only 2 parameters
-to_check = [b, c]
-
-ode = ODESystem(eqs, t, name = :GoodwinOsc)
+to_check = [ode.b, ode.c]
 
-global_id = assess_identifiability(ode, measured_quantities = measured_quantities,
-    funcs_to_check = to_check, p = 0.9)
+global_id = assess_identifiability(ode, funcs_to_check = to_check, p = 0.9)
 ```
 
 Both parameters `b, c` are globally identifiable with probability `0.9` in this case.
diff --git a/docs/src/tutorials/programmatically_generating.md b/docs/src/tutorials/programmatically_generating.md
new file mode 100644
index 0000000000..ca11243c84
--- /dev/null
+++ b/docs/src/tutorials/programmatically_generating.md
@@ -0,0 +1,106 @@
+# [Programmatically Generating and Scripting ODESystems](@id programmatically)
+
+In the following tutorial we will discuss how to programmatically generate `ODESystem`s.
+This is for cases where one is writing functions that generating `ODESystem`s, for example
+if implementing a reader which parses some file format to generate an `ODESystem` (for example,
+SBML), or for writing functions that transform an `ODESystem` (for example, if you write a
+function that log-transforms a variable in an `ODESystem`).
+
+## The Representation of a ModelingToolkit System
+
+ModelingToolkit is built on [Symbolics.jl](https://symbolics.juliasymbolics.org/dev/),
+a symbolic Computer Algebra System (CAS) developed in Julia. As such, all CAS functionality
+is available on ModelingToolkit systems, such as symbolic differentiation, Groebner basis
+calculations, and whatever else you can think of. Under the hood, all ModelingToolkit
+variables and expressions are Symbolics.jl variables and expressions. Thus when scripting
+a ModelingToolkit system, one simply needs to generate Symbolics.jl variables and equations
+as demonstrated in the Symbolics.jl documentation. This looks like:
+
+```@example scripting
+using Symbolics
+@variables t x(t) y(t) # Define variables
+D = Differential(t) # Define a differential operator
+eqs = [D(x) ~ y
+    D(y) ~ x] # Define an array of equations
+```
+
+## The Non-DSL (non-`@mtkmodel`) Way of Defining an ODESystem
+
+Using `@mtkmodel` is the preferred way of defining ODEs with MTK. However, let us
+look at how we can define the same system without `@mtkmodel`. This is useful for
+defining PDESystem etc.
+
+```@example scripting
+using ModelingToolkit
+
+@variables t x(t)   # independent and dependent variables
+@parameters τ       # parameters
+@constants h = 1    # constants
+D = Differential(t) # define an operator for the differentiation w.r.t. time
+eqs = [D(x) ~ (h - x) / τ] # create an array of equations
+
+# your first ODE, consisting of a single equation, indicated by ~
+@named fol_model = ODESystem(eqs, t)
+
+# Perform the standard transformations and mark the model complete
+# Note: Complete models cannot be subsystems of other models!
+fol_model = complete(structural_simplify(fol_model))
+```
+
+As you can see, generating an ODESystem is as simple as creating the array of equations
+and passing it to the `ODESystem` constructor.
+
+## Understanding the Difference Between the Julia Variable and the Symbolic Variable
+
+In the most basic usage of ModelingToolkit and Symbolics, the name of the Julia variable
+and the symbolic variable are the same. For example, when we do:
+
+```@example scripting
+@variables a
+```
+
+the name of the symbolic variable is `a` and same with the Julia variable. However, we can
+de-couple these by setting `a` to a new symbolic variable, for example:
+
+```@example scripting
+b = only(@variables(a))
+```
+
+Now the Julia variable `b` refers to the variable named `a`. However, the downside of this current
+approach is that it requires that the user writing the script knows the name `a` that they want to
+place to the variable. But what if for example we needed to get the variable's name from a file?
+
+To do this, one can interpolate a symbol into the `@variables` macro using `$`. For example:
+
+```@example scripting
+a = :c
+b = only(@variables($a))
+```
+
+In this example, `@variables($a)` created a variable named `c`, and set this variable to `b`.
+
+Variables are not the only thing with names. For example, when you build a system, it knows its name
+that name is used in the namespacing. In the standard usage, again the Julia variable and the
+symbolic name are made the same via:
+
+```@example scripting
+@named fol_model = ODESystem(eqs, t)
+```
+
+However, one can decouple these two properties by noting that `@named` is simply shorthand for the
+following:
+
+```@example scripting
+fol_model = ODESystem(eqs, t; name = :fol_model)
+```
+
+Thus if we had read a name from a file and wish to populate an `ODESystem` with said name, we could do:
+
+```@example scripting
+namesym = :name_from_file
+fol_model = ODESystem(eqs, t; name = namesym)
+```
+
+## Warning About Mutation
+
+Be advsied that it's never a good idea to mutate an `ODESystem`, or any `AbstractSystem`.
diff --git a/ext/MTKBifurcationKitExt.jl b/ext/MTKBifurcationKitExt.jl
new file mode 100644
index 0000000000..9a966ba55c
--- /dev/null
+++ b/ext/MTKBifurcationKitExt.jl
@@ -0,0 +1,40 @@
+module MTKBifurcationKitExt
+
+println("BifurcationKit extension loaded")
+
+### Preparations ###
+
+# Imports
+using ModelingToolkit, Setfield
+import BifurcationKit
+
+### Creates BifurcationProblem Overloads ###
+
+# When input is a NonlinearSystem.
+function BifurcationKit.BifurcationProblem(nsys::NonlinearSystem, u0_bif, ps, bif_par, args...; plot_var=nothing, record_from_solution=BifurcationKit.record_sol_default, jac=true, kwargs...)
+    # Creates F and J functions.
+    ofun = NonlinearFunction(nsys; jac=jac)
+    F = ofun.f
+    J = jac ? ofun.jac : nothing
+
+    # Computes bifurcation parameter and plot var indexes.
+    bif_idx = findfirst(isequal(bif_par), parameters(nsys))
+    if !isnothing(plot_var) 
+        plot_idx = findfirst(isequal(plot_var), states(nsys))
+        record_from_solution = (x, p) -> x[plot_idx]
+    end
+
+    # Converts the input state guess.
+    u0_bif = ModelingToolkit.varmap_to_vars(u0_bif, states(nsys))
+    ps = ModelingToolkit.varmap_to_vars(ps, parameters(nsys))
+
+    return BifurcationKit.BifurcationProblem(F, u0_bif, ps, (@lens _[bif_idx]), args...; record_from_solution = record_from_solution, J = J, kwargs...)
+end
+
+# When input is a ODESystem.
+function BifurcationKit.BifurcationProblem(osys::ODESystem, args...; kwargs...)
+    nsys = NonlinearSystem([0 ~ eq.rhs for eq in equations(osys)], states(osys), parameters(osys); name=osys.name)
+    return BifurcationKit.BifurcationProblem(nsys, args...; kwargs...)
+end
+
+end # module
diff --git a/ext/MTKDeepDiffsExt.jl b/ext/MTKDeepDiffsExt.jl
index 92bc9ba2b9..1d361f96a3 100644
--- a/ext/MTKDeepDiffsExt.jl
+++ b/ext/MTKDeepDiffsExt.jl
@@ -4,7 +4,7 @@ using DeepDiffs, ModelingToolkit
 using ModelingToolkit.BipartiteGraphs: Label,
     BipartiteAdjacencyList, unassigned,
     HighlightInt
-using ModelingToolkit.SystemStructures: SystemStructure,
+using ModelingToolkit: SystemStructure,
     MatchedSystemStructure,
     SystemStructurePrintMatrix
 
diff --git a/src/ModelingToolkit.jl b/src/ModelingToolkit.jl
index d5c4172340..de63a67837 100644
--- a/src/ModelingToolkit.jl
+++ b/src/ModelingToolkit.jl
@@ -51,8 +51,8 @@ using PrecompileTools, Reexport
     using MLStyle
 
     using Reexport
+    using Symbolics
     using Symbolics: degree
-    @reexport using Symbolics
     using Symbolics: _parse_vars, value, @derivatives, get_variables,
         exprs_occur_in, solve_for, build_expr, unwrap, wrap,
         VariableSource, getname, variable, Connection, connect,
@@ -70,9 +70,10 @@ using PrecompileTools, Reexport
     import OrdinaryDiffEq
 
     import Graphs: SimpleDiGraph, add_edge!, incidence_matrix
-
-    @reexport using UnPack
 end
+
+@reexport using Symbolics
+@reexport using UnPack
 RuntimeGeneratedFunctions.init(@__MODULE__)
 
 export @derivatives
@@ -156,7 +157,6 @@ include("systems/dependency_graphs.jl")
 include("clock.jl")
 include("discretedomain.jl")
 include("systems/systemstructure.jl")
-using .SystemStructures
 include("systems/clock_inference.jl")
 include("systems/systems.jl")
 
@@ -172,6 +172,14 @@ for S in subtypes(ModelingToolkit.AbstractSystem)
     @eval convert_system(::Type{<:$S}, sys::$S) = sys
 end
 
+PrecompileTools.@compile_workload begin
+    using ModelingToolkit
+    @variables t x(t)
+    D = Differential(t)
+    @named sys = ODESystem([D(x) ~ -x])
+    prob = ODEProblem(structural_simplify(sys), [x => 30.0], (0, 100), [], jac = true)
+end
+
 export AbstractTimeDependentSystem,
     AbstractTimeIndependentSystem,
     AbstractMultivariateSystem
@@ -193,7 +201,7 @@ export JumpProblem, DiscreteProblem
 export NonlinearSystem, OptimizationSystem, ConstraintsSystem
 export alias_elimination, flatten
 export connect, domain_connect, @connector, Connection, Flow, Stream, instream
-export @component, @mtkmodel
+export @component, @mtkmodel, @mtkbuild
 export isinput, isoutput, getbounds, hasbounds, isdisturbance, istunable, getdist, hasdist,
     tunable_parameters, isirreducible, getdescription, hasdescription, isbinaryvar,
     isintegervar
diff --git a/src/structural_transformation/StructuralTransformations.jl b/src/structural_transformation/StructuralTransformations.jl
index 47b140b7d1..7289df4232 100644
--- a/src/structural_transformation/StructuralTransformations.jl
+++ b/src/structural_transformation/StructuralTransformations.jl
@@ -23,14 +23,18 @@ using ModelingToolkit: ODESystem, AbstractSystem, var_from_nested_derivative, Di
     IncrementalCycleTracker, add_edge_checked!, topological_sort,
     invalidate_cache!, Substitutions, get_or_construct_tearing_state,
     filter_kwargs, lower_varname, setio, SparseMatrixCLIL,
-    fast_substitute, get_fullvars, has_equations
+    fast_substitute, get_fullvars, has_equations, observed
 
 using ModelingToolkit.BipartiteGraphs
 import .BipartiteGraphs: invview, complete
 import ModelingToolkit: var_derivative!, var_derivative_graph!
 using Graphs
-using ModelingToolkit.SystemStructures
-using ModelingToolkit.SystemStructures: algeqs, EquationsView
+using ModelingToolkit: algeqs, EquationsView,
+    SystemStructure, TransformationState, TearingState, structural_simplify!,
+    isdiffvar, isdervar, isalgvar, isdiffeq, algeqs, is_only_discrete,
+    dervars_range, diffvars_range, algvars_range,
+    DiffGraph, complete!,
+    get_fullvars, system_subset
 
 using ModelingToolkit.DiffEqBase
 using ModelingToolkit.StaticArrays
diff --git a/src/systems/abstractsystem.jl b/src/systems/abstractsystem.jl
index 1e4531a33b..37f4a49e23 100644
--- a/src/systems/abstractsystem.jl
+++ b/src/systems/abstractsystem.jl
@@ -229,7 +229,9 @@ for prop in [:eqs
     :metadata
     :gui_metadata
     :discrete_subsystems
-    :unknown_states]
+    :unknown_states
+    :split_idxs
+    :parent]
     fname1 = Symbol(:get_, prop)
     fname2 = Symbol(:has_, prop)
     @eval begin
@@ -306,6 +308,9 @@ function Base.propertynames(sys::AbstractSystem; private = false)
 end
 
 function Base.getproperty(sys::AbstractSystem, name::Symbol; namespace = !iscomplete(sys))
+    if has_parent(sys) && (parent = get_parent(sys); parent !== nothing)
+        sys = parent
+    end
     wrap(getvar(sys, name; namespace = namespace))
 end
 function getvar(sys::AbstractSystem, name::Symbol; namespace = !iscomplete(sys))
@@ -1179,6 +1184,15 @@ macro component(expr)
     esc(component_post_processing(expr, false))
 end
 
+macro mtkbuild(expr)
+    named_expr = ModelingToolkit.named_expr(expr)
+    name = named_expr.args[1]
+    esc(quote
+        $named_expr
+        $name = $structural_simplify($name)
+    end)
+end
+
 """
 $(SIGNATURES)
 
@@ -1273,14 +1287,34 @@ See also [`linearize`](@ref) which provides a higher-level interface.
 function linearization_function(sys::AbstractSystem, inputs,
     outputs; simplify = false,
     initialize = true,
+    op = Dict(),
+    p = DiffEqBase.NullParameters(),
+    zero_dummy_der = false,
     kwargs...)
-    sys, diff_idxs, alge_idxs, input_idxs = io_preprocessing(sys, inputs, outputs; simplify,
+    ssys, diff_idxs, alge_idxs, input_idxs = io_preprocessing(sys, inputs, outputs;
+        simplify,
         kwargs...)
+    if zero_dummy_der
+        dummyder = setdiff(states(ssys), states(sys))
+        defs = Dict(x => 0.0 for x in dummyder)
+        @set! ssys.defaults = merge(defs, defaults(ssys))
+        op = merge(defs, op)
+    end
+    sys = ssys
+    x0 = merge(defaults(sys), op)
+    u0, p, _ = get_u0_p(sys, x0, p; use_union = false, tofloat = true)
+    p, split_idxs = split_parameters_by_type(p)
+    ps = parameters(sys)
+    if p isa Tuple
+        ps = Base.Fix1(getindex, ps).(split_idxs)
+        ps = (ps...,) #if p is Tuple, ps should be Tuple
+    end
+
     lin_fun = let diff_idxs = diff_idxs,
         alge_idxs = alge_idxs,
         input_idxs = input_idxs,
         sts = states(sys),
-        fun = ODEFunction{true, SciMLBase.FullSpecialize}(sys),
+        fun = ODEFunction{true, SciMLBase.FullSpecialize}(sys, states(sys), ps; p = p),
         h = build_explicit_observed_function(sys, outputs),
         chunk = ForwardDiff.Chunk(input_idxs)
 
@@ -1599,11 +1633,12 @@ function linearize(sys, inputs, outputs; op = Dict(), t = 0.0,
     allow_input_derivatives = false,
     zero_dummy_der = false,
     kwargs...)
-    lin_fun, ssys = linearization_function(sys, inputs, outputs; kwargs...)
-    if zero_dummy_der
-        dummyder = setdiff(states(ssys), states(sys))
-        op = merge(op, Dict(x => 0.0 for x in dummyder))
-    end
+    lin_fun, ssys = linearization_function(sys,
+        inputs,
+        outputs;
+        zero_dummy_der,
+        op,
+        kwargs...)
     linearize(ssys, lin_fun; op, t, allow_input_derivatives), ssys
 end
 
diff --git a/src/systems/diffeqs/abstractodesystem.jl b/src/systems/diffeqs/abstractodesystem.jl
index 8528dd3c12..f11d1283c6 100644
--- a/src/systems/diffeqs/abstractodesystem.jl
+++ b/src/systems/diffeqs/abstractodesystem.jl
@@ -354,6 +354,7 @@ function DiffEqBase.ODEFunction{iip, specialize}(sys::AbstractODESystem, dvs = s
     checkbounds = false,
     sparsity = false,
     analytic = nothing,
+    split_idxs = nothing,
     kwargs...) where {iip, specialize}
     f_gen = generate_function(sys, dvs, ps; expression = Val{eval_expression},
         expression_module = eval_module, checkbounds = checkbounds,
@@ -508,6 +509,7 @@ function DiffEqBase.ODEFunction{iip, specialize}(sys::AbstractODESystem, dvs = s
         nothing
     end
 
+    @set! sys.split_idxs = split_idxs
     ODEFunction{iip, specialize}(f;
         sys = sys,
         jac = _jac === nothing ? nothing : _jac,
@@ -765,7 +767,7 @@ Take dictionaries with initial conditions and parameters and convert them to num
 """
 function get_u0_p(sys,
     u0map,
-    parammap;
+    parammap = nothing;
     use_union = true,
     tofloat = true,
     symbolic_u0 = false)
@@ -773,7 +775,9 @@ function get_u0_p(sys,
     ps = parameters(sys)
 
     defs = defaults(sys)
-    defs = mergedefaults(defs, parammap, ps)
+    if parammap !== nothing
+        defs = mergedefaults(defs, parammap, ps)
+    end
     defs = mergedefaults(defs, u0map, dvs)
 
     if symbolic_u0
@@ -835,7 +839,7 @@ function process_DEProblem(constructor, sys::AbstractODESystem, u0map, parammap;
     f = constructor(sys, dvs, ps, u0; ddvs = ddvs, tgrad = tgrad, jac = jac,
         checkbounds = checkbounds, p = p,
         linenumbers = linenumbers, parallel = parallel, simplify = simplify,
-        sparse = sparse, eval_expression = eval_expression,
+        sparse = sparse, eval_expression = eval_expression, split_idxs,
         kwargs...)
     implicit_dae ? (f, du0, u0, p) : (f, u0, p)
 end
diff --git a/src/systems/diffeqs/odesystem.jl b/src/systems/diffeqs/odesystem.jl
index dcbecfcfb0..9259953dfc 100644
--- a/src/systems/diffeqs/odesystem.jl
+++ b/src/systems/diffeqs/odesystem.jl
@@ -139,6 +139,14 @@ struct ODESystem <: AbstractODESystem
     used for ODAEProblem.
     """
     unknown_states::Union{Nothing, Vector{Any}}
+    """
+    split_idxs: a vector of vectors of indices for the split parameters.
+    """
+    split_idxs::Union{Nothing, Vector{Vector{Int}}}
+    """
+    parent: the hierarchical parent system before simplification.
+    """
+    parent::Any
 
     function ODESystem(tag, deqs, iv, dvs, ps, tspan, var_to_name, ctrls, observed, tgrad,
         jac, ctrl_jac, Wfact, Wfact_t, name, systems, defaults,
@@ -146,8 +154,8 @@ struct ODESystem <: AbstractODESystem
         devents, metadata = nothing, gui_metadata = nothing,
         tearing_state = nothing,
         substitutions = nothing, complete = false,
-        discrete_subsystems = nothing, unknown_states = nothing;
-        checks::Union{Bool, Int} = true)
+        discrete_subsystems = nothing, unknown_states = nothing,
+        split_idxs = nothing, parent = nothing; checks::Union{Bool, Int} = true)
         if checks == true || (checks & CheckComponents) > 0
             check_variables(dvs, iv)
             check_parameters(ps, iv)
@@ -161,7 +169,7 @@ struct ODESystem <: AbstractODESystem
             ctrl_jac, Wfact, Wfact_t, name, systems, defaults, torn_matching,
             connector_type, preface, cevents, devents, metadata, gui_metadata,
             tearing_state, substitutions, complete, discrete_subsystems,
-            unknown_states)
+            unknown_states, split_idxs, parent)
     end
 end
 
diff --git a/src/systems/diffeqs/sdesystem.jl b/src/systems/diffeqs/sdesystem.jl
index 30d0e10666..8c42f7ea0d 100644
--- a/src/systems/diffeqs/sdesystem.jl
+++ b/src/systems/diffeqs/sdesystem.jl
@@ -115,13 +115,17 @@ struct SDESystem <: AbstractODESystem
     complete: if a model `sys` is complete, then `sys.x` no longer performs namespacing.
     """
     complete::Bool
+    """
+    parent: the hierarchical parent system before simplification.
+    """
+    parent::Any
 
     function SDESystem(tag, deqs, neqs, iv, dvs, ps, tspan, var_to_name, ctrls, observed,
         tgrad,
         jac,
         ctrl_jac, Wfact, Wfact_t, name, systems, defaults, connector_type,
         cevents, devents, metadata = nothing, gui_metadata = nothing,
-        complete = false;
+        complete = false, parent = nothing;
         checks::Union{Bool, Int} = true)
         if checks == true || (checks & CheckComponents) > 0
             check_variables(dvs, iv)
@@ -135,7 +139,7 @@ struct SDESystem <: AbstractODESystem
         new(tag, deqs, neqs, iv, dvs, ps, tspan, var_to_name, ctrls, observed, tgrad, jac,
             ctrl_jac,
             Wfact, Wfact_t, name, systems, defaults, connector_type, cevents, devents,
-            metadata, gui_metadata, complete)
+            metadata, gui_metadata, complete, parent)
     end
 end
 
diff --git a/src/systems/discrete_system/discrete_system.jl b/src/systems/discrete_system/discrete_system.jl
index c0307c586b..ecd0cae8c6 100644
--- a/src/systems/discrete_system/discrete_system.jl
+++ b/src/systems/discrete_system/discrete_system.jl
@@ -86,6 +86,10 @@ struct DiscreteSystem <: AbstractTimeDependentSystem
     complete: if a model `sys` is complete, then `sys.x` no longer performs namespacing.
     """
     complete::Bool
+    """
+    parent: the hierarchical parent system before simplification.
+    """
+    parent::Any
 
     function DiscreteSystem(tag, discreteEqs, iv, dvs, ps, tspan, var_to_name, ctrls,
         observed,
@@ -93,7 +97,7 @@ struct DiscreteSystem <: AbstractTimeDependentSystem
         systems, defaults, preface, connector_type,
         metadata = nothing, gui_metadata = nothing,
         tearing_state = nothing, substitutions = nothing,
-        complete = false; checks::Union{Bool, Int} = true)
+        complete = false, parent = nothing; checks::Union{Bool, Int} = true)
         if checks == true || (checks & CheckComponents) > 0
             check_variables(dvs, iv)
             check_parameters(ps, iv)
@@ -105,7 +109,7 @@ struct DiscreteSystem <: AbstractTimeDependentSystem
             systems,
             defaults,
             preface, connector_type, metadata, gui_metadata,
-            tearing_state, substitutions, complete)
+            tearing_state, substitutions, complete, parent)
     end
 end
 
diff --git a/src/systems/jumps/jumpsystem.jl b/src/systems/jumps/jumpsystem.jl
index 8e87d145c2..b37aad5d19 100644
--- a/src/systems/jumps/jumpsystem.jl
+++ b/src/systems/jumps/jumpsystem.jl
@@ -293,7 +293,7 @@ dprob = DiscreteProblem(js, u₀map, tspan, parammap)
 function DiffEqBase.DiscreteProblem(sys::JumpSystem, u0map, tspan::Union{Tuple, Nothing},
     parammap = DiffEqBase.NullParameters();
     checkbounds = false,
-    use_union = false,
+    use_union = true,
     kwargs...)
     dvs = states(sys)
     ps = parameters(sys)
diff --git a/src/systems/model_parsing.jl b/src/systems/model_parsing.jl
index 6a95a6e9e8..299793c236 100644
--- a/src/systems/model_parsing.jl
+++ b/src/systems/model_parsing.jl
@@ -237,7 +237,7 @@ function parse_model!(exprs, comps, ext, eqs, icon, vs, ps, sps,
     if mname == Symbol("@components")
         parse_components!(exprs, comps, dict, body, kwargs)
     elseif mname == Symbol("@extend")
-        parse_extend!(exprs, ext, dict, body, kwargs)
+        parse_extend!(exprs, ext, dict, mod, body, kwargs)
     elseif mname == Symbol("@variables")
         parse_variables!(exprs, vs, dict, mod, body, :variables, kwargs)
     elseif mname == Symbol("@parameters")
@@ -372,7 +372,29 @@ function extend_args!(a, b, dict, expr, kwargs, varexpr, has_param = false)
     end
 end
 
-function parse_extend!(exprs, ext, dict, body, kwargs)
+const EMPTY_DICT = Dict()
+const EMPTY_VoVoSYMBOL = Vector{Symbol}[]
+
+function Base.names(model::Model)
+    vars = keys(get(model.structure, :variables, EMPTY_DICT))
+    vars = union(vars, keys(get(model.structure, :parameters, EMPTY_DICT)))
+    vars = union(vars,
+        map(first, get(model.structure, :components, EMPTY_VoVoSYMBOL)))
+    collect(vars)
+end
+
+function _parse_extend!(ext, a, b, dict, expr, kwargs, varexpr, vars)
+    extend_args!(a, b, dict, expr, kwargs, varexpr)
+    ext[] = a
+    push!(b.args, Expr(:kw, :name, Meta.quot(a)))
+    push!(expr.args, :($a = $b))
+
+    dict[:extend] = [Symbol.(vars.args), a, b.args[1]]
+
+    push!(expr.args, :(@unpack $vars = $a))
+end
+
+function parse_extend!(exprs, ext, dict, mod, body, kwargs)
     expr = Expr(:block)
     varexpr = Expr(:block)
     push!(exprs, varexpr)
@@ -380,28 +402,33 @@ function parse_extend!(exprs, ext, dict, body, kwargs)
     body = deepcopy(body)
     MLStyle.@match body begin
         Expr(:(=), a, b) => begin
-            vars = nothing
             if Meta.isexpr(b, :(=))
                 vars = a
                 if !Meta.isexpr(vars, :tuple)
                     error("`@extend` destructuring only takes an tuple as LHS. Got $body")
                 end
                 a, b = b.args
-                extend_args!(a, b, dict, expr, kwargs, varexpr)
-                vars, a, b
+                _parse_extend!(ext, a, b, dict, expr, kwargs, varexpr, vars)
+            else
+                error("When explicitly destructing in `@extend` please use the syntax: `@extend a, b = oneport = OnePort()`.")
             end
-            ext[] = a
-            push!(b.args, Expr(:kw, :name, Meta.quot(a)))
-            push!(expr.args, :($a = $b))
-
-            dict[:extend] = [Symbol.(vars.args), a, b.args[1]]
-
-            if vars !== nothing
-                push!(expr.args, :(@unpack $vars = $a))
+        end
+        Expr(:call, a′, _...) => begin
+            a = Symbol(Symbol("#mtkmodel"), :__anonymous__, a′)
+            b = body
+            if (model = getproperty(mod, b.args[1])) isa Model
+                vars = Expr(:tuple)
+                append!(vars.args, names(model))
+                _parse_extend!(ext, a, b, dict, expr, kwargs, varexpr, vars)
+            else
+                error("Cannot infer the exact `Model` that `@extend $(body)` refers." *
+                      " Please specify the names that it brings into scope by:" *
+                      " `@extend a, b = oneport = OnePort()`.")
             end
         end
         _ => error("`@extend` only takes an assignment expression. Got $body")
     end
+    return nothing
 end
 
 function parse_variable_arg!(expr, vs, dict, mod, arg, varclass, kwargs)
diff --git a/src/systems/nonlinear/nonlinearsystem.jl b/src/systems/nonlinear/nonlinearsystem.jl
index d33a198710..192089b7da 100644
--- a/src/systems/nonlinear/nonlinearsystem.jl
+++ b/src/systems/nonlinear/nonlinearsystem.jl
@@ -75,18 +75,23 @@ struct NonlinearSystem <: AbstractTimeIndependentSystem
     complete: if a model `sys` is complete, then `sys.x` no longer performs namespacing.
     """
     complete::Bool
+    """
+    parent: the hierarchical parent system before simplification.
+    """
+    parent::Any
 
     function NonlinearSystem(tag, eqs, states, ps, var_to_name, observed, jac, name,
         systems,
         defaults, connector_type, metadata = nothing,
         gui_metadata = nothing,
         tearing_state = nothing, substitutions = nothing,
-        complete = false; checks::Union{Bool, Int} = true)
+        complete = false, parent = nothing; checks::Union{Bool, Int} = true)
         if checks == true || (checks & CheckUnits) > 0
             all_dimensionless([states; ps]) || check_units(eqs)
         end
         new(tag, eqs, states, ps, var_to_name, observed, jac, name, systems, defaults,
-            connector_type, metadata, gui_metadata, tearing_state, substitutions, complete)
+            connector_type, metadata, gui_metadata, tearing_state, substitutions, complete,
+            parent)
     end
 end
 
diff --git a/src/systems/optimization/optimizationsystem.jl b/src/systems/optimization/optimizationsystem.jl
index 2d935aa030..7a9b10cd49 100644
--- a/src/systems/optimization/optimizationsystem.jl
+++ b/src/systems/optimization/optimizationsystem.jl
@@ -56,10 +56,14 @@ struct OptimizationSystem <: AbstractOptimizationSystem
     complete: if a model `sys` is complete, then `sys.x` no longer performs namespacing.
     """
     complete::Bool
+    """
+    parent: the hierarchical parent system before simplification.
+    """
+    parent::Any
 
     function OptimizationSystem(tag, op, states, ps, var_to_name, observed,
         constraints, name, systems, defaults, metadata = nothing,
-        gui_metadata = nothing, complete = false;
+        gui_metadata = nothing, complete = false, parent = nothing;
         checks::Union{Bool, Int} = true)
         if checks == true || (checks & CheckUnits) > 0
             unwrap(op) isa Symbolic && check_units(op)
@@ -67,7 +71,8 @@ struct OptimizationSystem <: AbstractOptimizationSystem
             all_dimensionless([states; ps]) || check_units(constraints)
         end
         new(tag, op, states, ps, var_to_name, observed,
-            constraints, name, systems, defaults, metadata, gui_metadata, complete)
+            constraints, name, systems, defaults, metadata, gui_metadata, complete,
+            parent)
     end
 end
 
diff --git a/src/systems/systems.jl b/src/systems/systems.jl
index 7aaaabb917..bdf8f31808 100644
--- a/src/systems/systems.jl
+++ b/src/systems/systems.jl
@@ -17,6 +17,23 @@ This will convert all `inputs` to parameters and allow them to be unconnected, i
 simplification will allow models where `n_states = n_equations - n_inputs`.
 """
 function structural_simplify(sys::AbstractSystem, io = nothing; simplify = false,
+    kwargs...)
+    newsys′ = __structural_simplify(sys, io; simplify, kwargs...)
+    if newsys′ isa Tuple
+        @assert length(newsys′) == 2
+        newsys = newsys′[1]
+    else
+        newsys = newsys′
+    end
+    @set! newsys.parent = complete(sys)
+    newsys = complete(newsys)
+    if newsys′ isa Tuple
+        return newsys, newsys′[2]
+    else
+        return newsys
+    end
+end
+function __structural_simplify(sys::AbstractSystem, io = nothing; simplify = false,
     kwargs...)
     sys = expand_connections(sys)
     sys isa DiscreteSystem && return sys
diff --git a/src/systems/systemstructure.jl b/src/systems/systemstructure.jl
index 4e330619d5..18e50afc9b 100644
--- a/src/systems/systemstructure.jl
+++ b/src/systems/systemstructure.jl
@@ -1,7 +1,5 @@
-module SystemStructures
-
 using DataStructures
-using Symbolics: linear_expansion, unwrap
+using Symbolics: linear_expansion, unwrap, Connection
 using SymbolicUtils: istree, operation, arguments, Symbolic
 using SymbolicUtils: quick_cancel, similarterm
 using ..ModelingToolkit
@@ -10,7 +8,7 @@ import ..ModelingToolkit: isdiffeq, var_from_nested_derivative, vars!, flatten,
     isparameter, isconstant,
     independent_variables, SparseMatrixCLIL, AbstractSystem,
     equations, isirreducible, input_timedomain, TimeDomain,
-    VariableType, getvariabletype, has_equations
+    VariableType, getvariabletype, has_equations, ODESystem
 using ..BipartiteGraphs
 import ..BipartiteGraphs: invview, complete
 using Graphs
@@ -27,8 +25,7 @@ function quick_cancel_expr(expr)
 end
 
 export SystemStructure, TransformationState, TearingState, structural_simplify!
-export initialize_system_structure, find_linear_equations
-export isdiffvar, isdervar, isalgvar, isdiffeq, isalgeq, algeqs, is_only_discrete
+export isdiffvar, isdervar, isalgvar, isdiffeq, algeqs, is_only_discrete
 export dervars_range, diffvars_range, algvars_range
 export DiffGraph, complete!
 export get_fullvars, system_subset
@@ -620,5 +617,3 @@ function _structural_simplify!(state::TearingState, io; simplify = false,
     @set! sys.observed = ModelingToolkit.topsort_equations(observed(sys), fullstates)
     ModelingToolkit.invalidate_cache!(sys), input_idxs
 end
-
-end # module
diff --git a/src/utils.jl b/src/utils.jl
index abea65ab21..ab17bd3d8d 100644
--- a/src/utils.jl
+++ b/src/utils.jl
@@ -657,46 +657,40 @@ function promote_to_concrete(vs; tofloat = true, use_union = true)
     end
     T = eltype(vs)
     if Base.isconcretetype(T) && (!tofloat || T === float(T)) # nothing to do
-        vs
+        return vs
     else
         sym_vs = filter(x -> SymbolicUtils.issym(x) || SymbolicUtils.istree(x), vs)
         isempty(sym_vs) || throw_missingvars_in_sys(sym_vs)
-        C = typeof(first(vs))
-        I = Int8
-        has_int = false
-        has_array = false
-        has_bool = false
-        array_T = nothing
+
+        C = nothing
         for v in vs
-            if v isa AbstractArray
-                has_array = true
-                array_T = typeof(v)
+            E = typeof(v)
+            if E <: Number
+                if tofloat
+                    E = float(E)
+                end
             end
-            E = eltype(v)
-            C = promote_type(C, E)
-            if E <: Integer
-                has_int = true
-                I = promote_type(I, E)
+            if C === nothing
+                C = E
             end
-            if E <: Bool
-                has_bool = true
+            if use_union
+                C = Union{C, E}
+            else
+                @assert C==E "`promote_to_concrete` can't make type $E uniform with $C"
+                C = E
             end
         end
-        if tofloat && !has_array
-            C = float(C)
-        elseif has_array || (use_union && has_int && C !== I)
-            if has_array
-                C = Union{C, array_T}
-            end
-            if has_int
-                C = Union{C, I}
-            end
-            if has_bool
-                C = Union{C, Bool}
+
+        y = similar(vs, C)
+        for i in eachindex(vs)
+            if (vs[i] isa Number) & tofloat
+                y[i] = float(vs[i]) #needed because copyto! can't convert Int to Float automatically
+            else
+                y[i] = vs[i]
             end
-            return copyto!(similar(vs, C), vs)
         end
-        convert.(C, vs)
+
+        return y
     end
 end
 
diff --git a/src/variables.jl b/src/variables.jl
index 4cffd3fdeb..854680998e 100644
--- a/src/variables.jl
+++ b/src/variables.jl
@@ -145,15 +145,23 @@ function SciMLBase.process_p_u0_symbolic(prob::Union{SciMLBase.AbstractDEProblem
                                 " Please use `remake` with the `u0` keyword argument as a vector of values, paying attention to state order."))
     end
 
-    # assemble defaults
-    defs = defaults(prob.f.sys)
-    defs = mergedefaults(defs, prob.p, parameters(prob.f.sys))
-    defs = mergedefaults(defs, p, parameters(prob.f.sys))
-    defs = mergedefaults(defs, prob.u0, states(prob.f.sys))
-    defs = mergedefaults(defs, u0, states(prob.f.sys))
-
-    u0 = varmap_to_vars(u0, states(prob.f.sys); defaults = defs, tofloat = true)
-    p = varmap_to_vars(p, parameters(prob.f.sys); defaults = defs)
+    sys = prob.f.sys
+    defs = defaults(sys)
+    ps = parameters(sys)
+    if has_split_idxs(sys) && (split_idxs = get_split_idxs(sys)) !== nothing
+        for (i, idxs) in enumerate(split_idxs)
+            defs = mergedefaults(defs, prob.p[i], ps[idxs])
+        end
+    else
+        # assemble defaults
+        defs = defaults(sys)
+        defs = mergedefaults(defs, prob.p, ps)
+    end
+    defs = mergedefaults(defs, p, ps)
+    sts = states(sys)
+    defs = mergedefaults(defs, prob.u0, sts)
+    defs = mergedefaults(defs, u0, sts)
+    u0, p, defs = get_u0_p(sys, defs)
 
     return p, u0
 end
diff --git a/test/clock.jl b/test/clock.jl
index 42074b65ff..74946d21d7 100644
--- a/test/clock.jl
+++ b/test/clock.jl
@@ -62,13 +62,12 @@ By inference:
 =>   Shift(x, 0, dt) := (Shift(x, -1, dt) + dt) / (1 - dt) # Discrete system
 =#
 
-using ModelingToolkit.SystemStructures
 ci, varmap = infer_clocks(sys)
 eqmap = ci.eq_domain
 tss, inputs = ModelingToolkit.split_system(deepcopy(ci))
-sss, = SystemStructures._structural_simplify!(deepcopy(tss[1]), (inputs[1], ()))
+sss, = ModelingToolkit._structural_simplify!(deepcopy(tss[1]), (inputs[1], ()))
 @test equations(sss) == [D(x) ~ u - x]
-sss, = SystemStructures._structural_simplify!(deepcopy(tss[2]), (inputs[2], ()))
+sss, = ModelingToolkit._structural_simplify!(deepcopy(tss[2]), (inputs[2], ()))
 @test isempty(equations(sss))
 @test observed(sss) == [yd ~ Sample(t, dt)(y); r ~ 1.0; ud ~ kp * (r - yd)]
 
diff --git a/test/extensions/bifurcationkit.jl b/test/extensions/bifurcationkit.jl
new file mode 100644
index 0000000000..c9ae8890ad
--- /dev/null
+++ b/test/extensions/bifurcationkit.jl
@@ -0,0 +1,29 @@
+using BifurcationKit, ModelingToolkit, Test
+
+# Checks pitchfork diagram and that there are the correct number of branches (a main one and two children)
+let 
+    @variables t x(t) y(t)
+    @parameters μ α
+    eqs = [0 ~ μ*x - x^3 + α*y,
+            0 ~ -y]
+    @named nsys = NonlinearSystem(eqs, [x, y], [μ, α])
+
+    bif_par = μ
+    p_start = [μ => -1.0, α => 1.0]
+    u0_guess = [x => 1.0, y => 1.0]
+    plot_var = x;
+
+    using BifurcationKit
+    bprob = BifurcationProblem(nsys, u0_guess, p_start, bif_par; plot_var=plot_var, jac=false)
+
+    p_span = (-4.0, 6.0)
+    opt_newton = NewtonPar(tol = 1e-9, max_iterations = 20)
+    opts_br = ContinuationPar(dsmin = 0.001, dsmax = 0.05, ds = 0.01,
+        max_steps = 100, nev = 2, newton_options = opt_newton,
+        p_min = p_span[1], p_max = p_span[2],
+        detect_bifurcation = 3, n_inversion = 4, tol_bisection_eigenvalue = 1e-8, dsmin_bisection = 1e-9);
+
+    bf = bifurcationdiagram(bprob, PALC(), 2, (args...) -> opts_br; bothside=true)
+
+    @test length(bf.child) == 2
+end
\ No newline at end of file
diff --git a/test/input_output_handling.jl b/test/input_output_handling.jl
index 12f68f9f04..249a94e9d0 100644
--- a/test/input_output_handling.jl
+++ b/test/input_output_handling.jl
@@ -124,8 +124,8 @@ using ModelingToolkitStandardLibrary.Mechanical.Rotational
 t = ModelingToolkitStandardLibrary.Mechanical.Rotational.t
 @named inertia1 = Inertia(; J = 1)
 @named inertia2 = Inertia(; J = 1)
-@named spring = Spring(; c = 10)
-@named damper = Damper(; d = 3)
+@named spring = Rotational.Spring(; c = 10)
+@named damper = Rotational.Damper(; d = 3)
 @named torque = Torque(; use_support = false)
 @variables y(t) = 0
 eqs = [connect(torque.flange, inertia1.flange_a)
diff --git a/test/latexify/10.tex b/test/latexify/10.tex
index 9a8335bd34..09ec74ab72 100644
--- a/test/latexify/10.tex
+++ b/test/latexify/10.tex
@@ -1,5 +1,5 @@
 \begin{align}
-\frac{\mathrm{d} x\left( t \right)}{\mathrm{d}t} =& \frac{\sigma \left(  - x\left( t \right) + y\left( t \right) \right) \frac{\mathrm{d}}{\mathrm{d}t} \left(  - y\left( t \right) + x\left( t \right) \right)}{\frac{\mathrm{d} z\left( t \right)}{\mathrm{d}t}} \\
-0 =&  - y\left( t \right) + \frac{1}{10} \sigma \left( \rho - z\left( t \right) \right) x\left( t \right) \\
-\frac{\mathrm{d} z\left( t \right)}{\mathrm{d}t} =& \left( y\left( t \right) \right)^{\frac{2}{3}} x\left( t \right) - \beta z\left( t \right)
+\frac{\mathrm{d} x\left( t \right)}{\mathrm{d}t} =& \frac{\left(  - x\left( t \right) + y\left( t \right) \right) \frac{\mathrm{d}}{\mathrm{d}t} \left( x\left( t \right) - y\left( t \right) \right) \sigma}{\frac{\mathrm{d} z\left( t \right)}{\mathrm{d}t}} \\
+0 =&  - y\left( t \right) + \frac{1}{10} x\left( t \right) \left(  - z\left( t \right) + \rho \right) \sigma \\
+\frac{\mathrm{d} z\left( t \right)}{\mathrm{d}t} =& \left( y\left( t \right) \right)^{\frac{2}{3}} x\left( t \right) - z\left( t \right) \beta
 \end{align}
diff --git a/test/latexify/20.tex b/test/latexify/20.tex
index e32ea81eb7..9de213aafd 100644
--- a/test/latexify/20.tex
+++ b/test/latexify/20.tex
@@ -1,5 +1,5 @@
 \begin{align}
-\frac{\mathrm{d}}{\mathrm{d}t} u(t)_1 =& \left(  - u(t)_1 + u(t)_2 \right) p_3 \\
-0 =&  - u(t)_2 + \frac{1}{10} \left(  - u(t)_1 + p_1 \right) p_2 p_3 u(t)_1 \\
+\frac{\mathrm{d}}{\mathrm{d}t} u(t)_1 =& p_3 \left(  - u(t)_1 + u(t)_2 \right) \\
+0 =&  - u(t)_2 + \frac{1}{10} \left( p_1 - u(t)_1 \right) p_2 p_3 u(t)_1 \\
 \frac{\mathrm{d}}{\mathrm{d}t} u(t)_3 =& u(t)_2^{\frac{2}{3}} u(t)_1 - p_3 u(t)_3
 \end{align}
diff --git a/test/latexify/30.tex b/test/latexify/30.tex
index bc46f41e9b..b83feeba72 100644
--- a/test/latexify/30.tex
+++ b/test/latexify/30.tex
@@ -1,5 +1,5 @@
 \begin{align}
-\frac{\mathrm{d}}{\mathrm{d}t} u(t)_1 =& \left(  - u(t)_1 + u(t)_2 \right) p_3 \\
-\frac{\mathrm{d}}{\mathrm{d}t} u(t)_2 =&  - u(t)_2 + \frac{1}{10} \left(  - u(t)_1 + p_1 \right) p_2 p_3 u(t)_1 \\
+\frac{\mathrm{d}}{\mathrm{d}t} u(t)_1 =& p_3 \left(  - u(t)_1 + u(t)_2 \right) \\
+\frac{\mathrm{d}}{\mathrm{d}t} u(t)_2 =&  - u(t)_2 + \frac{1}{10} \left( p_1 - u(t)_1 \right) p_2 p_3 u(t)_1 \\
 \frac{\mathrm{d}}{\mathrm{d}t} u(t)_3 =& u(t)_2^{\frac{2}{3}} u(t)_1 - p_3 u(t)_3
 \end{align}
diff --git a/test/model_parsing.jl b/test/model_parsing.jl
index 3f380df100..701e785928 100644
--- a/test/model_parsing.jl
+++ b/test/model_parsing.jl
@@ -1,83 +1,82 @@
 using ModelingToolkit, Test
 using ModelingToolkit: get_gui_metadata,
-    VariableDescription, getdefault, RegularConnector, get_ps
+    VariableDescription, getdefault, RegularConnector, get_ps, getname
 using URIs: URI
 using Distributions
 using Unitful
 
 ENV["MTK_ICONS_DIR"] = "$(@__DIR__)/icons"
 
-@testset "Comprehensive Test of Parsing Models (with an RC Circuit)" begin
-    @connector RealInput begin
-        u(t), [input = true, unit = u"V"]
+@connector RealInput begin
+    u(t), [input = true, unit = u"V"]
+end
+@connector RealOutput begin
+    u(t), [output = true, unit = u"V"]
+end
+@mtkmodel Constant begin
+    @components begin
+        output = RealOutput()
     end
-    @connector RealOutput begin
-        u(t), [output = true, unit = u"V"]
+    @parameters begin
+        k, [description = "Constant output value of block"]
     end
-    @mtkmodel Constant begin
-        @components begin
-            output = RealOutput()
-        end
-        @parameters begin
-            k, [description = "Constant output value of block"]
-        end
-        @equations begin
-            output.u ~ k
-        end
+    @equations begin
+        output.u ~ k
     end
+end
 
-    @variables t [unit = u"s"]
-    D = Differential(t)
+@variables t [unit = u"s"]
+D = Differential(t)
 
-    @connector Pin begin
-        v(t), [unit = u"V"]                    # Potential at the pin [V]
-        i(t), [connect = Flow, unit = u"A"]    # Current flowing into the pin [A]
-        @icon "pin.png"
-    end
+@connector Pin begin
+    v(t), [unit = u"V"]                    # Potential at the pin [V]
+    i(t), [connect = Flow, unit = u"A"]    # Current flowing into the pin [A]
+    @icon "pin.png"
+end
 
-    @named p = Pin(; v = π)
-    @test getdefault(p.v) == π
-    @test Pin.isconnector == true
+@named p = Pin(; v = π)
+@test getdefault(p.v) == π
+@test Pin.isconnector == true
 
-    @mtkmodel OnePort begin
-        @components begin
-            p = Pin()
-            n = Pin()
-        end
-        @variables begin
-            v(t), [unit = u"V"]
-            i(t), [unit = u"A"]
-        end
-        @icon "oneport.png"
-        @equations begin
-            v ~ p.v - n.v
-            0 ~ p.i + n.i
-            i ~ p.i
-        end
+@mtkmodel OnePort begin
+    @components begin
+        p = Pin()
+        n = Pin()
+    end
+    @variables begin
+        v(t), [unit = u"V"]
+        i(t), [unit = u"A"]
     end
+    @icon "oneport.png"
+    @equations begin
+        v ~ p.v - n.v
+        0 ~ p.i + n.i
+        i ~ p.i
+    end
+end
 
-    @test OnePort.isconnector == false
+@test OnePort.isconnector == false
 
-    @mtkmodel Ground begin
-        @components begin
-            g = Pin()
-        end
-        @icon begin
-            read(abspath(ENV["MTK_ICONS_DIR"], "ground.svg"), String)
-        end
-        @equations begin
-            g.v ~ 0
-        end
+@mtkmodel Ground begin
+    @components begin
+        g = Pin()
+    end
+    @icon begin
+        read(abspath(ENV["MTK_ICONS_DIR"], "ground.svg"), String)
     end
+    @equations begin
+        g.v ~ 0
+    end
+end
 
-    resistor_log = "$(@__DIR__)/logo/resistor.svg"
-    @mtkmodel Resistor begin
-        @extend v, i = oneport = OnePort()
-        @parameters begin
-            R, [unit = u"Ω"]
-        end
-        @icon begin
-            """<?xml version="1.0" encoding="UTF-8"?>
+resistor_log = "$(@__DIR__)/logo/resistor.svg"
+@mtkmodel Resistor begin
+    @extend v, i = oneport = OnePort()
+    @parameters begin
+        R, [unit = u"Ω"]
+    end
+    @icon begin
+        """<?xml version="1.0" encoding="UTF-8"?>
 <svg xmlns="http://www.w3.org/2000/svg" width="80" height="30">
   <path d="M10 15
 l15 0
@@ -91,84 +90,87 @@ l2.5 -5
 l15 0" stroke="black" stroke-width="1" stroke-linejoin="bevel" fill="none"></path>
 </svg>
 """
-        end
-        @equations begin
-            v ~ i * R
-        end
     end
+    @equations begin
+        v ~ i * R
+    end
+end
 
-    @mtkmodel Capacitor begin
-        @parameters begin
-            C, [unit = u"F"]
-        end
-        @extend v, i = oneport = OnePort(; v = 0.0)
-        @icon "https://upload.wikimedia.org/wikipedia/commons/7/78/Capacitor_symbol.svg"
-        @equations begin
-            D(v) ~ i / C
-        end
+@mtkmodel Capacitor begin
+    @parameters begin
+        C, [unit = u"F"]
+    end
+    @extend OnePort(; v = 0.0)
+    @icon "https://upload.wikimedia.org/wikipedia/commons/7/78/Capacitor_symbol.svg"
+    @equations begin
+        D(v) ~ i / C
     end
+end
 
-    @named capacitor = Capacitor(C = 10, v = 10.0)
-    @test getdefault(capacitor.v) == 10.0
+@named capacitor = Capacitor(C = 10, v = 10.0)
+@test getdefault(capacitor.v) == 10.0
 
-    @mtkmodel Voltage begin
-        @extend v, i = oneport = OnePort()
-        @components begin
-            V = RealInput()
-        end
-        @equations begin
-            v ~ V.u
-        end
+@mtkmodel Voltage begin
+    @extend v, i = oneport = OnePort()
+    @components begin
+        V = RealInput()
     end
+    @equations begin
+        v ~ V.u
+    end
+end
 
-    @mtkmodel RC begin
-        @structural_parameters begin
-            R_val = 10
-            C_val = 10
-            k_val = 10
-        end
-        @components begin
-            resistor = Resistor(; R = R_val)
-            capacitor = Capacitor(; C = C_val)
-            source = Voltage()
-            constant = Constant(; k = k_val)
-            ground = Ground()
-        end
-
-        @equations begin
-            connect(constant.output, source.V)
-            connect(source.p, resistor.p)
-            connect(resistor.n, capacitor.p)
-            connect(capacitor.n, source.n, ground.g)
-        end
+@mtkmodel RC begin
+    @structural_parameters begin
+        R_val = 10
+        C_val = 10
+        k_val = 10
+    end
+    @components begin
+        resistor = Resistor(; R = R_val)
+        capacitor = Capacitor(; C = C_val)
+        source = Voltage()
+        constant = Constant(; k = k_val)
+        ground = Ground()
     end
 
-    C_val = 20
-    R_val = 20
-    res__R = 100
-    @named rc = RC(; C_val, R_val, resistor.R = res__R)
-    # Test that `resistor.R` overrides `R_val` in the argument.
-    @test getdefault(rc.resistor.R) == res__R != R_val
-    # Test that `C_val` passed via argument is set as default of C.
-    @test getdefault(rc.capacitor.C) == C_val
-    # Test that `k`'s default value is unchanged.
-    @test getdefault(rc.constant.k) == RC.structure[:kwargs][:k_val]
-    @test getdefault(rc.capacitor.v) == 0.0
-
-    @test get_gui_metadata(rc.resistor).layout == Resistor.structure[:icon] ==
-          read(joinpath(ENV["MTK_ICONS_DIR"], "resistor.svg"), String)
-    @test get_gui_metadata(rc.ground).layout ==
-          read(abspath(ENV["MTK_ICONS_DIR"], "ground.svg"), String)
-    @test get_gui_metadata(rc.capacitor).layout ==
-          URI("https://upload.wikimedia.org/wikipedia/commons/7/78/Capacitor_symbol.svg")
-    @test OnePort.structure[:icon] ==
-          URI("file:///" * abspath(ENV["MTK_ICONS_DIR"], "oneport.png"))
-    @test ModelingToolkit.get_gui_metadata(rc.resistor.p).layout == Pin.structure[:icon] ==
-          URI("file:///" * abspath(ENV["MTK_ICONS_DIR"], "pin.png"))
-
-    @test length(equations(structural_simplify(rc))) == 1
+    @equations begin
+        connect(constant.output, source.V)
+        connect(source.p, resistor.p)
+        connect(resistor.n, capacitor.p)
+        connect(capacitor.n, source.n, ground.g)
+    end
 end
 
+C_val = 20
+R_val = 20
+res__R = 100
+@mtkbuild rc = RC(; C_val, R_val, resistor.R = res__R)
+resistor = getproperty(rc, :resistor; namespace = false)
+@test getname(rc.resistor) === getname(resistor)
+@test getname(rc.resistor.R) === getname(resistor.R)
+@test getname(rc.resistor.v) === getname(resistor.v)
+# Test that `resistor.R` overrides `R_val` in the argument.
+@test getdefault(rc.resistor.R) == res__R != R_val
+# Test that `C_val` passed via argument is set as default of C.
+@test getdefault(rc.capacitor.C) == C_val
+# Test that `k`'s default value is unchanged.
+@test getdefault(rc.constant.k) == RC.structure[:kwargs][:k_val]
+@test getdefault(rc.capacitor.v) == 0.0
+
+@test get_gui_metadata(rc.resistor).layout == Resistor.structure[:icon] ==
+      read(joinpath(ENV["MTK_ICONS_DIR"], "resistor.svg"), String)
+@test get_gui_metadata(rc.ground).layout ==
+      read(abspath(ENV["MTK_ICONS_DIR"], "ground.svg"), String)
+@test get_gui_metadata(rc.capacitor).layout ==
+      URI("https://upload.wikimedia.org/wikipedia/commons/7/78/Capacitor_symbol.svg")
+@test OnePort.structure[:icon] ==
+      URI("file:///" * abspath(ENV["MTK_ICONS_DIR"], "oneport.png"))
+@test ModelingToolkit.get_gui_metadata(rc.resistor.p).layout == Pin.structure[:icon] ==
+      URI("file:///" * abspath(ENV["MTK_ICONS_DIR"], "pin.png"))
+
+@test length(equations(rc)) == 1
+
 @testset "Parameters and Structural parameters in various modes" begin
     @mtkmodel MockModel begin
         @parameters begin
diff --git a/test/odesystem.jl b/test/odesystem.jl
index 56d436d76c..39bcda58de 100644
--- a/test/odesystem.jl
+++ b/test/odesystem.jl
@@ -23,7 +23,8 @@ eqs = [D(x) ~ σ * (y - x),
 ModelingToolkit.toexpr.(eqs)[1]
 @named de = ODESystem(eqs; defaults = Dict(x => 1))
 subed = substitute(de, [σ => k])
-@test isequal(sort(parameters(subed), by = string), [k, β, ρ])
+ssort(eqs) = sort(eqs, by = string)
+@test isequal(ssort(parameters(subed)), [k, β, ρ])
 @test isequal(equations(subed),
     [D(x) ~ k * (y - x)
         D(y) ~ (ρ - z) * x - y
@@ -241,7 +242,7 @@ p2 = (k₁ => 0.04,
     k₃ => 1e4)
 tspan = (0.0, 100000.0)
 prob1 = ODEProblem(sys, u0, tspan, p)
-@test prob1.f.sys === sys
+@test prob1.f.sys == sys
 prob12 = ODEProblem(sys, u0, tspan, [0.04, 3e7, 1e4])
 prob13 = ODEProblem(sys, u0, tspan, (0.04, 3e7, 1e4))
 prob14 = ODEProblem(sys, u0, tspan, p2)
@@ -348,14 +349,14 @@ eqs = [0 ~ x + z
     D(accumulation_y) ~ y
     D(accumulation_z) ~ z
     D(x) ~ y]
-@test sort(equations(asys), by = string) == eqs
+@test ssort(equations(asys)) == ssort(eqs)
 @variables ac(t)
 asys = add_accumulations(sys, [ac => (x + y)^2])
 eqs = [0 ~ x + z
     0 ~ x - y
     D(ac) ~ (x + y)^2
     D(x) ~ y]
-@test sort(equations(asys), by = string) == eqs
+@test ssort(equations(asys)) == ssort(eqs)
 
 sys2 = ode_order_lowering(sys)
 M = ModelingToolkit.calculate_massmatrix(sys2)
diff --git a/test/runtests.jl b/test/runtests.jl
index e6387e2cd9..8e1c124afa 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -55,5 +55,8 @@ end
 @safetestset "OptimizationSystem Test" include("optimizationsystem.jl")
 @safetestset "FuncAffect Test" include("funcaffect.jl")
 @safetestset "Constants Test" include("constants.jl")
+@safetestset "BifurcationKit Extension Test" include("extensions/bifurcationkit.jl")
 # Reference tests go Last
-@safetestset "Latexify recipes Test" include("latexify.jl")
+if VERSION >= v"1.9"
+    @safetestset "Latexify recipes Test" include("latexify.jl")
+end
diff --git a/test/split_parameters.jl b/test/split_parameters.jl
index ef8f434dca..d37b9e2c13 100644
--- a/test/split_parameters.jl
+++ b/test/split_parameters.jl
@@ -2,6 +2,29 @@ using ModelingToolkit, Test
 using ModelingToolkitStandardLibrary.Blocks
 using OrdinaryDiffEq
 
+x = [1, 2.0, false, [1, 2, 3], Parameter(1.0)]
+
+y = ModelingToolkit.promote_to_concrete(x)
+@test eltype(y) == Union{Float64, Parameter{Float64}, Vector{Int64}}
+
+y = ModelingToolkit.promote_to_concrete(x; tofloat = false)
+@test eltype(y) == Union{Bool, Float64, Int64, Parameter{Float64}, Vector{Int64}}
+
+x = [1, 2.0, false, [1, 2, 3]]
+y = ModelingToolkit.promote_to_concrete(x)
+@test eltype(y) == Union{Float64, Vector{Int64}}
+
+x = Any[1, 2.0, false]
+y = ModelingToolkit.promote_to_concrete(x; tofloat = false)
+@test eltype(y) == Union{Bool, Float64, Int64}
+
+y = ModelingToolkit.promote_to_concrete(x; use_union = false)
+@test eltype(y) == Float64
+
+x = Float16[1.0, 2.0, 3.0]
+y = ModelingToolkit.promote_to_concrete(x)
+@test eltype(y) == Float16
+
 # ------------------------ Mixed Single Values and Vector
 
 dt = 4e-4
@@ -46,12 +69,25 @@ eqs = [y ~ src.output.u
 @named sys = ODESystem(eqs, t, vars, []; systems = [int, src])
 s = complete(sys)
 sys = structural_simplify(sys)
-prob = ODEProblem(sys, [], (0.0, t_end), [s.src.data => x])
+prob = ODEProblem(sys, [], (0.0, t_end), [s.src.data => x]; tofloat = false)
 @test prob.p isa Tuple{Vector{Float64}, Vector{Int}, Vector{Vector{Float64}}}
 sol = solve(prob, ImplicitEuler());
 @test sol.retcode == ReturnCode.Success
 @test sol[y][end] == x[end]
 
+#TODO: remake becomes more complicated now, how to improve?
+defs = ModelingToolkit.defaults(sys)
+defs[s.src.data] = 2x
+p′ = ModelingToolkit.varmap_to_vars(defs, parameters(sys); tofloat = false)
+p′, = ModelingToolkit.split_parameters_by_type(p′) #NOTE: we need to ensure this is called now before calling remake()
+prob′ = remake(prob; p = p′)
+sol = solve(prob′, ImplicitEuler());
+@test sol.retcode == ReturnCode.Success
+@test sol[y][end] == 2x[end]
+
+prob′′ = remake(prob; p = [s.src.data => x])
+@test_broken prob′′.p isa Tuple
+
 # ------------------------ Mixed Type Converted to float (default behavior)
 
 vars = @variables y(t)=1 dy(t)=0 ddy(t)=0
@@ -77,3 +113,74 @@ prob = ODEProblem(sys, [], tspan, []; tofloat = false)
 @test prob.p isa Tuple{Vector{Float64}, Vector{Int64}}
 sol = solve(prob, ImplicitEuler());
 @test sol.retcode == ReturnCode.Success
+
+# ------------------------- Bug
+using ModelingToolkit, LinearAlgebra
+using ModelingToolkitStandardLibrary.Mechanical.Rotational
+using ModelingToolkitStandardLibrary.Blocks
+using ModelingToolkitStandardLibrary.Blocks: t
+using ModelingToolkit: connect
+
+"A wrapper function to make symbolic indexing easier"
+function wr(sys)
+    ODESystem(Equation[], ModelingToolkit.get_iv(sys), systems = [sys], name = :a_wrapper)
+end
+indexof(sym, syms) = findfirst(isequal(sym), syms)
+
+# Parameters
+m1 = 1.0
+m2 = 1.0
+k = 10.0 # Spring stiffness
+c = 3.0  # Damping coefficient
+
+@named inertia1 = Inertia(; J = m1)
+@named inertia2 = Inertia(; J = m2)
+@named spring = Spring(; c = k)
+@named damper = Damper(; d = c)
+@named torque = Torque(use_support = false)
+
+function SystemModel(u = nothing; name = :model)
+    eqs = [connect(torque.flange, inertia1.flange_a)
+        connect(inertia1.flange_b, spring.flange_a, damper.flange_a)
+        connect(inertia2.flange_a, spring.flange_b, damper.flange_b)]
+    if u !== nothing
+        push!(eqs, connect(torque.tau, u.output))
+        return @named model = ODESystem(eqs,
+            t;
+            systems = [torque, inertia1, inertia2, spring, damper, u])
+    end
+    ODESystem(eqs, t; systems = [torque, inertia1, inertia2, spring, damper], name)
+end
+
+model = SystemModel() # Model with load disturbance
+@named d = Step(start_time = 1.0, duration = 10.0, offset = 0.0, height = 1.0) # Disturbance
+model_outputs = [model.inertia1.w, model.inertia2.w, model.inertia1.phi, model.inertia2.phi] # This is the state realization we want to control
+inputs = [model.torque.tau.u]
+matrices, ssys = ModelingToolkit.linearize(wr(model), inputs, model_outputs)
+
+# Design state-feedback gain using LQR
+# Define cost matrices
+x_costs = [model.inertia1.w => 1.0
+    model.inertia2.w => 1.0
+    model.inertia1.phi => 1.0
+    model.inertia2.phi => 1.0]
+L = randn(1, 4) # Post-multiply by `C` to get the correct input to the controller
+
+# This old definition of MatrixGain will work because the parameter space does not include K (an Array term)
+# @component function MatrixGainAlt(K::AbstractArray; name)
+#     nout, nin = size(K, 1), size(K, 2)
+#     @named input = RealInput(; nin = nin)
+#     @named output = RealOutput(; nout = nout)
+#     eqs = [output.u[i] ~ sum(K[i, j] * input.u[j] for j in 1:nin) for i in 1:nout]
+#     compose(ODESystem(eqs, t, [], []; name = name), [input, output])
+# end
+
+@named state_feedback = MatrixGain(K = -L) # Build negative feedback into the feedback matrix
+@named add = Add(; k1 = 1.0, k2 = 1.0) # To add the control signal and the disturbance
+
+connections = [[state_feedback.input.u[i] ~ model_outputs[i] for i in 1:4]
+    connect(d.output, :d, add.input1)
+    connect(add.input2, state_feedback.output)
+    connect(add.output, :u, model.torque.tau)]
+@named closed_loop = ODESystem(connections, t, systems = [model, state_feedback, add, d])
+S = get_sensitivity(closed_loop, :u)
diff --git a/test/stream_connectors.jl b/test/stream_connectors.jl
index 9502b55be9..e4365876cf 100644
--- a/test/stream_connectors.jl
+++ b/test/stream_connectors.jl
@@ -136,29 +136,30 @@ eqns = [domain_connect(fluid, n1m1.port_a)
 
 @test_nowarn structural_simplify(n1m1Test)
 @unpack source, port_a = n1m1
-@test sort(equations(expand_connections(n1m1)), by = string) == [0 ~ port_a.m_flow
+ssort(eqs) = sort(eqs, by = string)
+@test ssort(equations(expand_connections(n1m1))) == ssort([0 ~ port_a.m_flow
     0 ~ source.port1.m_flow - port_a.m_flow
     source.port1.P ~ port_a.P
     source.port1.P ~ source.P
     source.port1.h_outflow ~ port_a.h_outflow
-    source.port1.h_outflow ~ source.h]
+    source.port1.h_outflow ~ source.h])
 @unpack port_a, port_b = pipe
-@test sort(equations(expand_connections(pipe)), by = string) ==
-      [0 ~ -port_a.m_flow - port_b.m_flow
+@test ssort(equations(expand_connections(pipe))) ==
+      ssort([0 ~ -port_a.m_flow - port_b.m_flow
     0 ~ port_a.m_flow
     0 ~ port_b.m_flow
     port_a.P ~ port_b.P
     port_a.h_outflow ~ instream(port_b.h_outflow)
-    port_b.h_outflow ~ instream(port_a.h_outflow)]
-@test sort(equations(expand_connections(sys)), by = string) ==
-      [0 ~ n1m1.port_a.m_flow + pipe.port_a.m_flow
+    port_b.h_outflow ~ instream(port_a.h_outflow)])
+@test ssort(equations(expand_connections(sys))) ==
+      ssort([0 ~ n1m1.port_a.m_flow + pipe.port_a.m_flow
     0 ~ pipe.port_b.m_flow + sink.port.m_flow
     n1m1.port_a.P ~ pipe.port_a.P
-    pipe.port_b.P ~ sink.port.P]
-@test sort(equations(expand_connections(n1m1Test)), by = string) ==
-      [0 ~ -pipe.port_a.m_flow - pipe.port_b.m_flow
-    0 ~ n1m1.port_a.m_flow + pipe.port_a.m_flow
+    pipe.port_b.P ~ sink.port.P])
+@test ssort(equations(expand_connections(n1m1Test))) ==
+      ssort([0 ~ -pipe.port_a.m_flow - pipe.port_b.m_flow
     0 ~ n1m1.source.port1.m_flow - n1m1.port_a.m_flow
+    0 ~ n1m1.port_a.m_flow + pipe.port_a.m_flow
     0 ~ pipe.port_b.m_flow + sink.port.m_flow
     fluid.m_flow ~ 0
     n1m1.port_a.P ~ pipe.port_a.P
@@ -172,7 +173,7 @@ eqns = [domain_connect(fluid, n1m1.port_a)
     pipe.port_b.h_outflow ~ n1m1.port_a.h_outflow
     sink.port.P ~ sink.P
     sink.port.h_outflow ~ sink.h_in
-    sink.port.m_flow ~ -sink.m_flow_in]
+    sink.port.m_flow ~ -sink.m_flow_in])
 
 # N1M2 model and test code.
 function N1M2(; name,
@@ -255,12 +256,12 @@ eqns = [connect(source.port, n2m2.port_a)
 @named sp2 = TwoPhaseFluidPort()
 @named sys = ODESystem([connect(sp1, sp2)], t)
 sys_exp = expand_connections(compose(sys, [sp1, sp2]))
-@test sort(equations(sys_exp), by = string) == [0 ~ -sp1.m_flow - sp2.m_flow
+@test ssort(equations(sys_exp)) == ssort([0 ~ -sp1.m_flow - sp2.m_flow
     0 ~ sp1.m_flow
     0 ~ sp2.m_flow
     sp1.P ~ sp2.P
     sp1.h_outflow ~ ModelingToolkit.instream(sp2.h_outflow)
-    sp2.h_outflow ~ ModelingToolkit.instream(sp1.h_outflow)]
+    sp2.h_outflow ~ ModelingToolkit.instream(sp1.h_outflow)])
 
 # array var
 @connector function VecPin(; name)
@@ -274,15 +275,15 @@ end
 
 @named simple = ODESystem([connect(vp1, vp2, vp3)], t)
 sys = expand_connections(compose(simple, [vp1, vp2, vp3]))
-@test sort(equations(sys), by = string) == sort([0 .~ collect(vp1.i)
-        0 .~ collect(vp2.i)
-        0 .~ collect(vp3.i)
-        vp1.v[1] ~ vp2.v[1]
-        vp1.v[2] ~ vp2.v[2]
-        vp1.v[1] ~ vp3.v[1]
-        vp1.v[2] ~ vp3.v[2]
-        0 ~ -vp1.i[1] - vp2.i[1] - vp3.i[1]
-        0 ~ -vp1.i[2] - vp2.i[2] - vp3.i[2]], by = string)
+@test ssort(equations(sys)) == ssort([0 .~ collect(vp1.i)
+    0 .~ collect(vp2.i)
+    0 .~ collect(vp3.i)
+    vp1.v[1] ~ vp2.v[1]
+    vp1.v[2] ~ vp2.v[2]
+    vp1.v[1] ~ vp3.v[1]
+    vp1.v[2] ~ vp3.v[2]
+    0 ~ -vp1.i[1] - vp2.i[1] - vp3.i[1]
+    0 ~ -vp1.i[2] - vp2.i[2] - vp3.i[2]])
 
 @connector function VectorHeatPort(; name, N = 100, T0 = 0.0, Q0 = 0.0)
     @variables (T(t))[1:N]=T0 (Q(t))[1:N]=Q0 [connect = Flow]