diff --git a/NEWS.md b/NEWS.md
index 94b8c507ca..1a1c6c4126 100644
--- a/NEWS.md
+++ b/NEWS.md
@@ -49,3 +49,4 @@
     `(single_parameter => expression_involving_other_parameters)` and a `Vector` of these can be passed to
     the `parameter_dependencies` keyword argument of `ODESystem`, `SDESystem` and `JumpSystem`. The dependent
     parameters are updated whenever other parameters are modified, e.g. in callbacks.
+  - Support for `IfElse.jl` has been dropped. `Base.ifelse` can be used instead.
diff --git a/Project.toml b/Project.toml
index b1f9d68908..aa0466e91d 100644
--- a/Project.toml
+++ b/Project.toml
@@ -23,7 +23,6 @@ FindFirstFunctions = "64ca27bc-2ba2-4a57-88aa-44e436879224"
 ForwardDiff = "f6369f11-7733-5829-9624-2563aa707210"
 FunctionWrappersWrappers = "77dc65aa-8811-40c2-897b-53d922fa7daf"
 Graphs = "86223c79-3864-5bf0-83f7-82e725a168b6"
-IfElse = "615f187c-cbe4-4ef1-ba3b-2fcf58d6d173"
 InteractiveUtils = "b77e0a4c-d291-57a0-90e8-8db25a27a240"
 JuliaFormatter = "98e50ef6-434e-11e9-1051-2b60c6c9e899"
 JumpProcesses = "ccbc3e58-028d-4f4c-8cd5-9ae44345cda5"
@@ -81,7 +80,6 @@ FindFirstFunctions = "1"
 ForwardDiff = "0.10.3"
 FunctionWrappersWrappers = "0.1"
 Graphs = "1.5.2"
-IfElse = "0.1"
 InteractiveUtils = "1"
 JuliaFormatter = "1.0.47"
 JumpProcesses = "9.1"
@@ -116,7 +114,6 @@ julia = "1.9"
 AmplNLWriter = "7c4d4715-977e-5154-bfe0-e096adeac482"
 BenchmarkTools = "6e4b80f9-dd63-53aa-95a3-0cdb28fa8baf"
 ControlSystemsBase = "aaaaaaaa-a6ca-5380-bf3e-84a91bcd477e"
-ControlSystemsMTK = "687d7614-c7e5-45fc-bfc3-9ee385575c88"
 DeepDiffs = "ab62b9b5-e342-54a8-a765-a90f495de1a6"
 ForwardDiff = "f6369f11-7733-5829-9624-2563aa707210"
 Ipopt = "b6b21f68-93f8-5de0-b562-5493be1d77c9"
@@ -139,4 +136,4 @@ Sundials = "c3572dad-4567-51f8-b174-8c6c989267f4"
 Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
 
 [targets]
-test = ["AmplNLWriter", "BenchmarkTools", "ControlSystemsBase", "ControlSystemsMTK", "NonlinearSolve", "ForwardDiff", "Ipopt", "Ipopt_jll", "ModelingToolkitStandardLibrary", "Optimization", "OptimizationOptimJL", "OptimizationMOI", "Random", "ReferenceTests", "SafeTestsets", "StableRNGs", "Statistics", "SteadyStateDiffEq", "Test", "StochasticDiffEq", "Sundials", "StochasticDelayDiffEq", "Pkg"]
+test = ["AmplNLWriter", "BenchmarkTools", "ControlSystemsBase", "NonlinearSolve", "ForwardDiff", "Ipopt", "Ipopt_jll", "ModelingToolkitStandardLibrary", "Optimization", "OptimizationOptimJL", "OptimizationMOI", "Random", "ReferenceTests", "SafeTestsets", "StableRNGs", "Statistics", "SteadyStateDiffEq", "Test", "StochasticDiffEq", "Sundials", "StochasticDelayDiffEq", "Pkg"]
diff --git a/docs/Project.toml b/docs/Project.toml
index 49fba20527..8aaabc6054 100644
--- a/docs/Project.toml
+++ b/docs/Project.toml
@@ -4,9 +4,9 @@ BifurcationKit = "0f109fa4-8a5d-4b75-95aa-f515264e7665"
 DifferentialEquations = "0c46a032-eb83-5123-abaf-570d42b7fbaa"
 Distributions = "31c24e10-a181-5473-b8eb-7969acd0382f"
 Documenter = "e30172f5-a6a5-5a46-863b-614d45cd2de4"
+DynamicQuantities = "06fc5a27-2a28-4c7c-a15d-362465fb6821"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 ModelingToolkit = "961ee093-0014-501f-94e3-6117800e7a78"
-ModelingToolkitDesigner = "23d639d0-9462-4d1e-84fe-d700424865b8"
 NonlinearSolve = "8913a72c-1f9b-4ce2-8d82-65094dcecaec"
 Optim = "429524aa-4258-5aef-a3af-852621145aeb"
 Optimization = "7f7a1694-90dd-40f0-9382-eb1efda571ba"
@@ -14,7 +14,7 @@ OptimizationOptimJL = "36348300-93cb-4f02-beb5-3c3902f8871e"
 OrdinaryDiffEq = "1dea7af3-3e70-54e6-95c3-0bf5283fa5ed"
 Plots = "91a5bcdd-55d7-5caf-9e0b-520d859cae80"
 StochasticDiffEq = "789caeaf-c7a9-5a7d-9973-96adeb23e2a0"
-StructuralIdentifiability = "220ca800-aa68-49bb-acd8-6037fa93a544"
+SymbolicIndexingInterface = "2efcf032-c050-4f8e-a9bb-153293bab1f5"
 SymbolicUtils = "d1185830-fcd6-423d-90d6-eec64667417b"
 Symbolics = "0c5d862f-8b57-4792-8d23-62f2024744c7"
 Unitful = "1986cc42-f94f-5a68-af5c-568840ba703d"
@@ -26,7 +26,7 @@ DifferentialEquations = "7.6"
 Distributions = "0.25"
 Documenter = "1"
 ModelingToolkit = "8.33, 9"
-ModelingToolkitDesigner = "1"
+DynamicQuantities = "^0.11.2"
 NonlinearSolve = "0.3, 1, 2, 3"
 Optim = "1.7"
 Optimization = "3.9"
@@ -34,7 +34,7 @@ OptimizationOptimJL = "0.1"
 OrdinaryDiffEq = "6.31"
 Plots = "1.36"
 StochasticDiffEq = "6"
-StructuralIdentifiability = "0.4, 0.5"
+SymbolicIndexingInterface = "0.3.1"
 SymbolicUtils = "1"
 Symbolics = "5"
 Unitful = "1.12"
diff --git a/docs/src/basics/Events.md b/docs/src/basics/Events.md
index 91df993f70..61d59863e2 100644
--- a/docs/src/basics/Events.md
+++ b/docs/src/basics/Events.md
@@ -144,12 +144,12 @@ ModelingToolkit therefore supports regular Julia functions as affects: instead
 of one or more equations, an affect is defined as a `tuple`:
 
 ```julia
-[x ~ 0] => (affect!, [v, x], [p, q], ctx)
+[x ~ 0] => (affect!, [v, x], [p, q], [discretes...], ctx)
 ```
 
 where, `affect!` is a Julia function with the signature: `affect!(integ, u, p, ctx)`; `[u,v]` and `[p,q]` are the symbolic unknowns (variables) and parameters
-that are accessed by `affect!`, respectively; and `ctx` is any context that is
-passed to `affect!` as the `ctx` argument.
+that are accessed by `affect!`, respectively; `discretes` are the parameters modified by `affect!`, if any;
+and `ctx` is any context that is passed to `affect!` as the `ctx` argument.
 
 `affect!` receives a [DifferentialEquations.jl
 integrator](https://docs.sciml.ai/DiffEqDocs/stable/basics/integrator/)
@@ -172,7 +172,7 @@ When accessing variables of a sub-system, it can be useful to rename them
 (alternatively, an affect function may be reused in different contexts):
 
 ```julia
-[x ~ 0] => (affect!, [resistor₊v => :v, x], [p, q => :p2], ctx)
+[x ~ 0] => (affect!, [resistor₊v => :v, x], [p, q => :p2], [], ctx)
 ```
 
 Here, the symbolic variable `resistor₊v` is passed as `v` while the symbolic
@@ -191,7 +191,7 @@ function bb_affect!(integ, u, p, ctx)
     integ.u[u.v] = -integ.u[u.v]
 end
 
-reflect = [x ~ 0] => (bb_affect!, [v], [], nothing)
+reflect = [x ~ 0] => (bb_affect!, [v], [], [], nothing)
 
 @mtkbuild bb_sys = ODESystem(bb_eqs, t, sts, par,
     continuous_events = reflect)
diff --git a/docs/src/basics/FAQ.md b/docs/src/basics/FAQ.md
index 19fcf68080..9d12b619fb 100644
--- a/docs/src/basics/FAQ.md
+++ b/docs/src/basics/FAQ.md
@@ -24,12 +24,11 @@ pnew = varmap_to_vars([β => 3.0, c => 10.0, γ => 2.0], parameters(sys))
 
 ## How do I handle `if` statements in my symbolic forms?
 
-For statements that are in the `if then else` form, use `IfElse.ifelse` from the
-[IfElse.jl](https://github.com/SciML/IfElse.jl) package to represent the code in a
-functional form. For handling direct `if` statements, you can use equivalent boolean
-mathematical expressions. For example, `if x > 0 ...` can be implemented as just
-`(x > 0) * `, where if `x <= 0` then the boolean will evaluate to `0` and thus the
-term will be excluded from the model.
+For statements that are in the `if then else` form, use `Base.ifelse` from the
+to represent the code in a functional form. For handling direct `if` statements,
+you can use equivalent boolean mathematical expressions. For example, `if x > 0 ...`
+can be implemented as just `(x > 0) * `, where if `x <= 0` then the boolean will
+evaluate to `0` and thus the term will be excluded from the model.
 
 ## ERROR: TypeError: non-boolean (Num) used in boolean context?
 
diff --git a/docs/src/basics/Validation.md b/docs/src/basics/Validation.md
index 346110f6d5..c9c662f13b 100644
--- a/docs/src/basics/Validation.md
+++ b/docs/src/basics/Validation.md
@@ -7,7 +7,7 @@ ModelingToolkit.jl provides extensive functionality for model validation and uni
 Units may be assigned with the following syntax.
 
 ```@example validation
-using ModelingToolkit, Unitful
+using ModelingToolkit, DynamicQuantities
 @variables t [unit = u"s"] x(t) [unit = u"m"] g(t) w(t) [unit = "Hz"]
 
 @variables(t, [unit = u"s"], x(t), [unit = u"m"], g(t), w(t), [unit = "Hz"])
@@ -45,7 +45,7 @@ ModelingToolkit.get_unit
 Example usage below. Note that `ModelingToolkit` does not force unit conversions to preferred units in the event of nonstandard combinations -- it merely checks that the equations are consistent.
 
 ```@example validation
-using ModelingToolkit, Unitful
+using ModelingToolkit, DynamicQuantities
 @parameters τ [unit = u"ms"]
 @variables t [unit = u"ms"] E(t) [unit = u"kJ"] P(t) [unit = u"MW"]
 D = Differential(t)
@@ -59,13 +59,17 @@ ModelingToolkit.validate(eqs[1])
 ```
 
 ```@example validation
-ModelingToolkit.get_unit(eqs[1].rhs)
+try
+    ModelingToolkit.get_unit(eqs[1].rhs)
+catch e
+    showerror(stdout, e)
+end
 ```
 
 An example of an inconsistent system: at present, `ModelingToolkit` requires that the units of all terms in an equation or sum to be equal-valued (`ModelingToolkit.equivalent(u1,u2)`), rather than simply dimensionally consistent. In the future, the validation stage may be upgraded to support the insertion of conversion factors into the equations.
 
 ```@example validation
-using ModelingToolkit, Unitful
+using ModelingToolkit, DynamicQuantities
 @parameters τ [unit = u"ms"]
 @variables t [unit = u"ms"] E(t) [unit = u"J"] P(t) [unit = u"MW"]
 D = Differential(t)
@@ -81,10 +85,9 @@ expect an object type, while two-parameter calls expect a function type as the f
 second argument.
 
 ```@example validation2
-using ModelingToolkit, Unitful
+using ModelingToolkit, DynamicQuantities
+using ModelingToolkit: t_nounits as t, D_nounits as D
 # Composite type parameter in registered function
-@variables t
-D = Differential(t)
 struct NewType
     f::Any
 end
@@ -101,16 +104,17 @@ end
 sts = @variables a(t)=0 [unit = u"cm"]
 ps = @parameters s=-1 [unit = u"cm"] c=c [unit = u"cm"]
 eqs = [D(a) ~ dummycomplex(c, s);]
-sys = ODESystem(eqs, t, [sts...;], [ps...;], name = :sys)
+sys = ODESystem(
+    eqs, t, [sts...;], [ps...;], name = :sys, checks = ~ModelingToolkit.CheckUnits)
 sys_simple = structural_simplify(sys)
 ```
 
-## `Unitful` Literals
+## `DynamicQuantities` Literals
 
 In order for a function to work correctly during both validation & execution, the function must be unit-agnostic. That is, no unitful literals may be used. Any unitful quantity must either be a `parameter` or `variable`. For example, these equations will not validate successfully.
 
 ```julia
-using ModelingToolkit, Unitful
+using ModelingToolkit, DynamicQuantities
 @variables t [unit = u"ms"] E(t) [unit = u"J"] P(t) [unit = u"MW"]
 D = Differential(t)
 eqs = [D(E) ~ P - E / 1u"ms"]
@@ -124,7 +128,7 @@ ModelingToolkit.validate(eqs) #Returns false while displaying a warning message
 Instead, they should be parameterized:
 
 ```@example validation3
-using ModelingToolkit, Unitful
+using ModelingToolkit, DynamicQuantities
 @parameters τ [unit = u"ms"]
 @variables t [unit = u"ms"] E(t) [unit = u"kJ"] P(t) [unit = u"MW"]
 D = Differential(t)
@@ -138,8 +142,8 @@ eqs = [D(E) ~ P - myfunc(E, τ)]
 ModelingToolkit.validate(eqs) #Returns true
 ```
 
-It is recommended *not* to circumvent unit validation by specializing user-defined functions on `Unitful` arguments vs. `Numbers`. This both fails to take advantage of `validate` for ensuring correctness, and may cause in errors in the
-future when `ModelingToolkit` is extended to support eliminating `Unitful` literals from functions.
+It is recommended *not* to circumvent unit validation by specializing user-defined functions on `DynamicQuantities` arguments vs. `Numbers`. This both fails to take advantage of `validate` for ensuring correctness, and may cause in errors in the
+future when `ModelingToolkit` is extended to support eliminating `DynamicQuantities` literals from functions.
 
 ## Other Restrictions
 
@@ -149,7 +153,7 @@ future when `ModelingToolkit` is extended to support eliminating `Unitful` liter
 
 If a system fails to validate due to unit issues, at least one warning message will appear, including a line number as well as the unit types and expressions that were in conflict. Some system constructors re-order equations before the unit checking can be done, in which case the equation numbers may be inaccurate. The printed expression that the problem resides in is always correctly shown.
 
-Symbolic exponents for unitful variables *are* supported (ex: `P^γ` in thermodynamics). However, this means that `ModelingToolkit` cannot reduce such expressions to `Unitful.Unitlike` subtypes at validation time because the exponent value is not available. In this case `ModelingToolkit.get_unit` is type-unstable, yielding a symbolic result, which can still be checked for symbolic equality with `ModelingToolkit.equivalent`.
+Symbolic exponents for unitful variables *are* supported (ex: `P^γ` in thermodynamics). However, this means that `ModelingToolkit` cannot reduce such expressions to `DynamicQuantities.Quantity` subtypes at validation time because the exponent value is not available. In this case `ModelingToolkit.get_unit` is type-unstable, yielding a symbolic result, which can still be checked for symbolic equality with `ModelingToolkit.equivalent`.
 
 ## Parameter & Initial Condition Values
 
diff --git a/docs/src/basics/Variable_metadata.md b/docs/src/basics/Variable_metadata.md
index a973832538..5e25f33373 100644
--- a/docs/src/basics/Variable_metadata.md
+++ b/docs/src/basics/Variable_metadata.md
@@ -142,7 +142,7 @@ In the example below, we define a system with tunable parameters and extract bou
 @variables x(t)=0 u(t)=0 [input = true] y(t)=0 [output = true]
 @parameters T [tunable = true, bounds = (0, Inf)]
 @parameters k [tunable = true, bounds = (0, Inf)]
-eqs = [Dₜ(x) ~ (-x + k * u) / T # A first-order system with time constant T and gain k
+eqs = [D(x) ~ (-x + k * u) / T # A first-order system with time constant T and gain k
        y ~ x]
 sys = ODESystem(eqs, t, name = :tunable_first_order)
 ```
@@ -159,6 +159,9 @@ lb, ub = getbounds(p) # operating on a vector, we get lower and upper bound vect
 b = getbounds(sys) # Operating on the system, we get a dict
 ```
 
+See also: [`ModelingToolkit.dump_variable_metadata`](@ref), [`ModelingToolkit.dump_parameters`](@ref),
+[`ModelingToolkit.dump_unknowns`](@ref).
+
 ## Index
 
 ```@index
@@ -172,3 +175,9 @@ Modules = [ModelingToolkit]
 Pages = ["variables.jl"]
 Private = false
 ```
+
+```@docs
+ModelingToolkit.dump_variable_metadata
+ModelingToolkit.dump_parameters
+ModelingToolkit.dump_unknowns
+```
diff --git a/docs/src/examples/parsing.md b/docs/src/examples/parsing.md
index 78ad824d6e..66a4e4d82f 100644
--- a/docs/src/examples/parsing.md
+++ b/docs/src/examples/parsing.md
@@ -23,10 +23,11 @@ From there, we can use ModelingToolkit to transform the symbolic equations into
 nonlinear solve:
 
 ```@example parsing
-using ModelingToolkit, NonlinearSolve
+using ModelingToolkit, SymbolicIndexingInterface, NonlinearSolve
 vars = union(ModelingToolkit.vars.(eqs)...)
 @mtkbuild ns = NonlinearSystem(eqs, vars, [])
 
-prob = NonlinearProblem(ns, [1.0, 1.0, 1.0])
+varmap = Dict(SymbolicIndexingInterface.getname.(vars) .=> vars)
+prob = NonlinearProblem(ns, [varmap[:x] => 1.0, varmap[:y] => 1.0, varmap[:z] => 1.0])
 sol = solve(prob, NewtonRaphson())
 ```
diff --git a/docs/src/internals.md b/docs/src/internals.md
index 79ded5c57f..a0d9ebfdf6 100644
--- a/docs/src/internals.md
+++ b/docs/src/internals.md
@@ -32,6 +32,6 @@ The call chain typically looks like this, with the function names in the case of
 
  1. Problem constructor ([`ODEProblem`](@ref))
  2. Build an `DEFunction` ([`process_DEProblem`](@ref) -> [`ODEFunction`](@ref)
- 3. Write actual executable code ([`generate_function`](@ref))
+ 3. Write actual executable code ([`generate_function`](@ref) or [`generate_custom_function`](@ref))
 
 Apart from [`generate_function`](@ref), which generates the dynamics function, `ODEFunction` also builds functions for observed equations (`build_explicit_observed_function`) and Jacobians (`generate_jacobian`) etc. These are all stored in the `ODEFunction`.
diff --git a/docs/src/tutorials/domain_connections.md b/docs/src/tutorials/domain_connections.md
index b381668544..d6dc2d8781 100644
--- a/docs/src/tutorials/domain_connections.md
+++ b/docs/src/tutorials/domain_connections.md
@@ -94,7 +94,7 @@ end
 nothing #hide
 ```
 
-When the system is defined we can generate a fluid component and connect it to the system.  Here `fluid` is connected to the `src.port`, but it could also be connected to `vol.port`, any connection in the network is fine.  Note: we can visualize the system using `ModelingToolkitDesigner.jl`, where a dashed line is used to show the `fluid` connection to represent a domain connection that is only transporting parameters and not unknown variables.
+When the system is defined we can generate a fluid component and connect it to the system.  Here `fluid` is connected to the `src.port`, but it could also be connected to `vol.port`, any connection in the network is fine.
 
 ```@example domain
 @component function System(; name)
@@ -115,20 +115,6 @@ end
 nothing #hide
 ```
 
-```@setup domain
-# code to generate diagrams...
-# using ModelingToolkitDesigner 
-# path = raw"C:\Work\Assets\ModelingToolkit.jl\domain_connections" 
-# design = ODESystemDesign(odesys, path); 
-
-# using CairoMakie
-# CairoMakie.set_theme!(Theme(;fontsize=12))
-# fig = ModelingToolkitDesigner.view(design, false)
-# save(joinpath(path, "odesys.svg"), fig; resolution=(300,300))
-```
-
-![odesys](https://github.com/SciML/ModelingToolkit.jl/assets/40798837/d19fbcf4-781c-4743-87b7-30bed348ff98)
-
 To see how the domain works, we can examine the set parameter values for each of the ports `src.port` and `vol.port`.  First we assemble the system using `structural_simplify()` and then check the default value of `vol.port.ρ`, whichs points to the setter value `fluid₊ρ`.  Likewise, `src.port.ρ`, will also point to the setter value `fluid₊ρ`.  Therefore, there is now only 1 defined density value `fluid₊ρ` which sets the density for the connected network.
 
 ```@repl domain
@@ -195,14 +181,6 @@ end
 nothing #hide
 ```
 
-```@setup domain
-# design = ODESystemDesign(actsys2, path);
-# fig = ModelingToolkitDesigner.view(design, false)
-# save(joinpath(path, "actsys2.svg"), fig; resolution=(500,300))
-```
-
-![actsys2](https://github.com/SciML/ModelingToolkit.jl/assets/40798837/8ed50035-f6ac-48cb-a585-1ef415154a02)
-
 After running `structural_simplify()` on `actsys2`, the defaults will show that `act.port_a.ρ` points to `fluid_a₊ρ` and `act.port_b.ρ` points to `fluid_b₊ρ`.  This is a special case, in most cases a hydraulic system will have only 1 fluid, however this simple system has 2 separate domain networks.  Therefore, we can connect a single fluid to both networks.  This does not interfere with the mathematical equations of the system, since no unknown variables are connected.
 
 ```@example domain
@@ -227,14 +205,6 @@ end
 nothing #hide
 ```
 
-```@setup domain
-# design = ODESystemDesign(actsys1, path);
-# fig = ModelingToolkitDesigner.view(design, false)
-# save(joinpath(path, "actsys1.svg"), fig; resolution=(500,300))
-```
-
-![actsys1](https://github.com/SciML/ModelingToolkit.jl/assets/40798837/054404eb-dbb7-4b85-95c0-c9503d0c4d00)
-
 ## Special Connection Cases (`domain_connect()`)
 
 In some cases a component will be defined with 2 connectors of the same domain, but they are not connected.  For example the `Restrictor` defined here gives equations to define the behavior of how the 2 connectors `port_a` and `port_b` are physically connected.
@@ -282,14 +252,6 @@ end
 nothing #hide
 ```
 
-```@setup domain
-# design = ODESystemDesign(ressys, path);
-# fig = ModelingToolkitDesigner.view(design, false)
-# save(joinpath(path, "ressys.svg"), fig; resolution=(500,300))
-```
-
-![ressys](https://github.com/SciML/ModelingToolkit.jl/assets/40798837/3740f0e2-7324-4c1f-af8b-eba02cfece81)
-
 When `structural_simplify()` is applied to this system it can be seen that the defaults are missing for `res.port_b` and `vol.port`.
 
 ```@repl domain
diff --git a/docs/src/tutorials/ode_modeling.md b/docs/src/tutorials/ode_modeling.md
index e749f42362..ef354071d9 100644
--- a/docs/src/tutorials/ode_modeling.md
+++ b/docs/src/tutorials/ode_modeling.md
@@ -338,7 +338,7 @@ While defining the model `UnitstepFOLFactory`, an initial guess of 0.0 is assign
 Additionally, these initial guesses can be modified while creating instances of `UnitstepFOLFactory` by passing arguments.
 
 ```@example ode2
-@named fol = UnitstepFOLFactory(; x = 0.1)
+@mtkbuild fol = UnitstepFOLFactory(; x = 0.1)
 sol = ODEProblem(fol, [], (0.0, 5.0), []) |> solve
 ```
 
diff --git a/docs/src/tutorials/optimization.md b/docs/src/tutorials/optimization.md
index 3aa4a25596..9f71b592d8 100644
--- a/docs/src/tutorials/optimization.md
+++ b/docs/src/tutorials/optimization.md
@@ -6,13 +6,13 @@ Let's solve the classical _Rosenbrock function_ in two dimensions.
 
 First, we need to make some imports.
 
-```@example rosenbrock_2d
+```julia
 using ModelingToolkit, Optimization, OptimizationOptimJL
 ```
 
 Now we can define our optimization problem.
 
-```@example rosenbrock_2d
+```julia
 @variables begin
     x, [bounds = (-2.0, 2.0)]
     y, [bounds = (-1.0, 3.0)]
@@ -44,7 +44,7 @@ Every optimization problem consists of a set of _optimization variables_. In thi
 
 Next, the actual `OptimizationProblem` can be created. At this stage, an initial guess `u0` for the optimization variables needs to be provided via map, using the symbols from before. Concrete values for the parameters of the system can also be provided or changed. However, if the parameters have default values assigned, they are used automatically.
 
-```@example rosenbrock_2d
+```julia
 u0 = [x => 1.0
       y => 2.0]
 p = [a => 1.0
@@ -56,7 +56,7 @@ solve(prob, GradientDescent())
 
 ## Rosenbrock Function with Constraints
 
-```@example rosenbrock_2d_cstr
+```julia
 using ModelingToolkit, Optimization, OptimizationOptimJL
 
 @variables begin
diff --git a/docs/src/tutorials/parameter_identifiability.md b/docs/src/tutorials/parameter_identifiability.md
index bcd027dcf0..54c3d0f42d 100644
--- a/docs/src/tutorials/parameter_identifiability.md
+++ b/docs/src/tutorials/parameter_identifiability.md
@@ -28,7 +28,7 @@ To define the ode system in Julia, we use `ModelingToolkit.jl`.
 
 We first define the parameters, variables, differential equations and the output equations.
 
-```@example SI
+```julia
 using StructuralIdentifiability, ModelingToolkit
 using ModelingToolkit: t_nounits as t, D_nounits as D
 
@@ -66,7 +66,7 @@ end
 
 After that, we are ready to check the system for local identifiability:
 
-```@example SI
+```julia
 # query local identifiability
 # we pass the ode-system
 local_id_all = assess_local_identifiability(de, p = 0.99)
@@ -76,7 +76,7 @@ We can see that all unknowns (except $x_7$) and all parameters are locally ident
 
 Let's try to check specific parameters and their combinations
 
-```@example SI
+```julia
 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)
 ```
@@ -103,7 +103,7 @@ 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.
 
-```@example SI2
+```julia
 using StructuralIdentifiability, ModelingToolkit
 using ModelingToolkit: t_nounits as t, D_nounits as D
 
@@ -144,7 +144,7 @@ We can see that only parameters `a, g` are unidentifiable, and everything else c
 
 Let us consider the same system but with two inputs, and we will find out identifiability with probability `0.9` for parameters `c` and `b`:
 
-```@example SI3
+```julia
 using StructuralIdentifiability, ModelingToolkit
 using ModelingToolkit: t_nounits as t, D_nounits as D
 
diff --git a/docs/src/tutorials/programmatically_generating.md b/docs/src/tutorials/programmatically_generating.md
index 079f50632d..8b8b03027a 100644
--- a/docs/src/tutorials/programmatically_generating.md
+++ b/docs/src/tutorials/programmatically_generating.md
@@ -41,7 +41,7 @@ using ModelingToolkit: t_nounits as t, D_nounits as D
 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)
+@named model = ODESystem(eqs, t)
 
 # Perform the standard transformations and mark the model complete
 # Note: Complete models cannot be subsystems of other models!
diff --git a/docs/src/tutorials/stochastic_diffeq.md b/docs/src/tutorials/stochastic_diffeq.md
index 91ec05ef01..fd03d51810 100644
--- a/docs/src/tutorials/stochastic_diffeq.md
+++ b/docs/src/tutorials/stochastic_diffeq.md
@@ -5,7 +5,7 @@ to the model: randomness. This is a
 [stochastic differential equation](https://en.wikipedia.org/wiki/Stochastic_differential_equation)
 which has a deterministic (drift) component and a stochastic (diffusion)
 component. Let's take the Lorenz equation from the first tutorial and extend
-it to have multiplicative noise.
+it to have multiplicative noise by creating `@brownian` variables in the equations.
 
 ```@example SDE
 using ModelingToolkit, StochasticDiffEq
@@ -14,16 +14,12 @@ using ModelingToolkit: t_nounits as t, D_nounits as D
 # Define some variables
 @parameters σ ρ β
 @variables x(t) y(t) z(t)
+@brownian a
+eqs = [D(x) ~ σ * (y - x) + 0.1a * x,
+    D(y) ~ x * (ρ - z) - y + 0.1a * y,
+    D(z) ~ x * y - β * z + 0.1a * z]
 
-eqs = [D(x) ~ σ * (y - x),
-    D(y) ~ x * (ρ - z) - y,
-    D(z) ~ x * y - β * z]
-
-noiseeqs = [0.1 * x,
-    0.1 * y,
-    0.1 * z]
-
-@mtkbuild de = SDESystem(eqs, noiseeqs, t, [x, y, z], [σ, ρ, β])
+@mtkbuild de = System(eqs, t)
 
 u0map = [
     x => 1.0,
@@ -38,5 +34,5 @@ parammap = [
 ]
 
 prob = SDEProblem(de, u0map, (0.0, 100.0), parammap)
-sol = solve(prob, SOSRI())
+sol = solve(prob, LambaEulerHeun())
 ```
diff --git a/src/ModelingToolkit.jl b/src/ModelingToolkit.jl
index 224ed22724..b70074f6fa 100644
--- a/src/ModelingToolkit.jl
+++ b/src/ModelingToolkit.jl
@@ -27,7 +27,6 @@ using PrecompileTools, Reexport
     using DocStringExtensions
     using Base: RefValue
     using Combinatorics
-    import IfElse
     import Distributions
     import FunctionWrappersWrappers
     using URIs: URI
@@ -222,8 +221,7 @@ export connect, domain_connect, @connector, Connection, Flow, Stream, instream
 export @component, @mtkmodel, @mtkbuild
 export isinput, isoutput, getbounds, hasbounds, getguess, hasguess, isdisturbance,
        istunable, getdist, hasdist,
-       tunable_parameters, isirreducible, getdescription, hasdescription, isbinaryvar,
-       isintegervar
+       tunable_parameters, isirreducible, getdescription, hasdescription
 export ode_order_lowering, dae_order_lowering, liouville_transform
 export PDESystem
 export Differential, expand_derivatives, @derivatives
@@ -234,7 +232,7 @@ export independent_variable, equations, controls,
        observed, full_equations
 export structural_simplify, expand_connections, linearize, linearization_function
 
-export calculate_jacobian, generate_jacobian, generate_function
+export calculate_jacobian, generate_jacobian, generate_function, generate_custom_function
 export calculate_control_jacobian, generate_control_jacobian
 export calculate_tgrad, generate_tgrad
 export calculate_gradient, generate_gradient
diff --git a/src/systems/abstractsystem.jl b/src/systems/abstractsystem.jl
index 7bb3c84818..f95cb1fe6f 100644
--- a/src/systems/abstractsystem.jl
+++ b/src/systems/abstractsystem.jl
@@ -82,7 +82,7 @@ function calculate_hessian end
 
 """
 ```julia
-generate_tgrad(sys::AbstractTimeDependentSystem, dvs = unknowns(sys), ps = parameters(sys),
+generate_tgrad(sys::AbstractTimeDependentSystem, dvs = unknowns(sys), ps = full_parameters(sys),
                expression = Val{true}; kwargs...)
 ```
 
@@ -93,7 +93,7 @@ function generate_tgrad end
 
 """
 ```julia
-generate_gradient(sys::AbstractSystem, dvs = unknowns(sys), ps = parameters(sys),
+generate_gradient(sys::AbstractSystem, dvs = unknowns(sys), ps = full_parameters(sys),
                   expression = Val{true}; kwargs...)
 ```
 
@@ -104,7 +104,7 @@ function generate_gradient end
 
 """
 ```julia
-generate_jacobian(sys::AbstractSystem, dvs = unknowns(sys), ps = parameters(sys),
+generate_jacobian(sys::AbstractSystem, dvs = unknowns(sys), ps = full_parameters(sys),
                   expression = Val{true}; sparse = false, kwargs...)
 ```
 
@@ -115,7 +115,7 @@ function generate_jacobian end
 
 """
 ```julia
-generate_factorized_W(sys::AbstractSystem, dvs = unknowns(sys), ps = parameters(sys),
+generate_factorized_W(sys::AbstractSystem, dvs = unknowns(sys), ps = full_parameters(sys),
                       expression = Val{true}; sparse = false, kwargs...)
 ```
 
@@ -126,7 +126,7 @@ function generate_factorized_W end
 
 """
 ```julia
-generate_hessian(sys::AbstractSystem, dvs = unknowns(sys), ps = parameters(sys),
+generate_hessian(sys::AbstractSystem, dvs = unknowns(sys), ps = full_parameters(sys),
                  expression = Val{true}; sparse = false, kwargs...)
 ```
 
@@ -137,7 +137,7 @@ function generate_hessian end
 
 """
 ```julia
-generate_function(sys::AbstractSystem, dvs = unknowns(sys), ps = parameters(sys),
+generate_function(sys::AbstractSystem, dvs = unknowns(sys), ps = full_parameters(sys),
                   expression = Val{true}; kwargs...)
 ```
 
@@ -145,6 +145,168 @@ Generate a function to evaluate the system's equations.
 """
 function generate_function end
 
+"""
+```julia
+generate_custom_function(sys::AbstractSystem, exprs, dvs = unknowns(sys),
+                         ps = full_parameters(sys); kwargs...)
+```
+
+Generate a function to evaluate `exprs`. `exprs` is a symbolic expression or
+array of symbolic expression involving symbolic variables in `sys`. The symbolic variables
+may be subsetted using `dvs` and `ps`. All `kwargs` except `postprocess_fbody` and `states`
+are passed to the internal [`build_function`](@ref) call. The returned function can be called
+as `f(u, p, t)` or `f(du, u, p, t)` for time-dependent systems and `f(u, p)` or `f(du, u, p)`
+for time-independent systems. If `split=true` (the default) was passed to [`complete`](@ref),
+[`structural_simplify`](@ref) or [`@mtkbuild`](@ref), `p` is expected to be an `MTKParameters`
+object.
+"""
+function generate_custom_function(sys::AbstractSystem, exprs, dvs = unknowns(sys),
+        ps = parameters(sys); wrap_code = nothing, kwargs...)
+    if !iscomplete(sys)
+        error("A completed system is required. Call `complete` or `structural_simplify` on the system.")
+    end
+    p = reorder_parameters(sys, ps)
+    isscalar = !(exprs isa AbstractArray)
+    if wrap_code === nothing
+        wrap_code = isscalar ? identity : (identity, identity)
+    end
+    pre, sol_states = get_substitutions_and_solved_unknowns(sys)
+
+    if is_time_dependent(sys)
+        return build_function(exprs,
+            dvs,
+            p...,
+            get_iv(sys);
+            kwargs...,
+            postprocess_fbody = pre,
+            states = sol_states,
+            wrap_code = wrap_code .∘ wrap_mtkparameters(sys, isscalar) .∘
+                        wrap_array_vars(sys, exprs; dvs)
+        )
+    else
+        return build_function(exprs,
+            dvs,
+            p...;
+            kwargs...,
+            postprocess_fbody = pre,
+            states = sol_states,
+            wrap_code = wrap_code .∘ wrap_mtkparameters(sys, isscalar) .∘
+                        wrap_array_vars(sys, exprs; dvs)
+        )
+    end
+end
+
+function wrap_array_vars(sys::AbstractSystem, exprs; dvs = unknowns(sys))
+    isscalar = !(exprs isa AbstractArray)
+    allvars = if isscalar
+        Set(get_variables(exprs))
+    else
+        union(get_variables.(exprs)...)
+    end
+    array_vars = Dict{Any, AbstractArray{Int}}()
+    for (j, x) in enumerate(dvs)
+        if istree(x) && operation(x) == getindex
+            arg = arguments(x)[1]
+            arg in allvars || continue
+            inds = get!(() -> Int[], array_vars, arg)
+            push!(inds, j)
+        end
+    end
+    for (k, inds) in array_vars
+        if inds == (inds′ = inds[1]:inds[end])
+            array_vars[k] = inds′
+        end
+    end
+    if isscalar
+        function (expr)
+            Func(
+                expr.args,
+                [],
+                Let(
+                    [k ← :(view($(expr.args[1].name), $v)) for (k, v) in array_vars],
+                    expr.body,
+                    false
+                )
+            )
+        end
+    else
+        function (expr)
+            Func(
+                expr.args,
+                [],
+                Let(
+                    [k ← :(view($(expr.args[1].name), $v)) for (k, v) in array_vars],
+                    expr.body,
+                    false
+                )
+            )
+        end,
+        function (expr)
+            Func(
+                expr.args,
+                [],
+                Let(
+                    [k ← :(view($(expr.args[2].name), $v)) for (k, v) in array_vars],
+                    expr.body,
+                    false
+                )
+            )
+        end
+    end
+end
+
+function wrap_mtkparameters(sys::AbstractSystem, isscalar::Bool)
+    if has_index_cache(sys) && get_index_cache(sys) !== nothing
+        offset = Int(is_time_dependent(sys))
+
+        if isscalar
+            function (expr)
+                p = gensym(:p)
+                Func(
+                    [
+                        expr.args[1],
+                        DestructuredArgs(
+                            [arg.name for arg in expr.args[2:(end - offset)]], p),
+                        (isone(offset) ? (expr.args[end],) : ())...
+                    ],
+                    [],
+                    Let(expr.args[2:(end - offset)], expr.body, false)
+                )
+            end
+        else
+            function (expr)
+                p = gensym(:p)
+                Func(
+                    [
+                        expr.args[1],
+                        DestructuredArgs(
+                            [arg.name for arg in expr.args[2:(end - offset)]], p),
+                        (isone(offset) ? (expr.args[end],) : ())...
+                    ],
+                    [],
+                    Let(expr.args[2:(end - offset)], expr.body, false)
+                )
+            end,
+            function (expr)
+                p = gensym(:p)
+                Func(
+                    [
+                        expr.args[1],
+                        expr.args[2],
+                        DestructuredArgs(
+                            [arg.name for arg in expr.args[3:(end - offset)]], p),
+                        (isone(offset) ? (expr.args[end],) : ())...
+                    ],
+                    [],
+                    Let(expr.args[3:(end - offset)], expr.body, false)
+                )
+            end
+        end
+    else
+        identity
+    end
+end
+
 mutable struct Substitutions
     subs::Vector{Equation}
     deps::Vector{Vector{Int}}
@@ -202,6 +364,9 @@ function SymbolicIndexingInterface.is_variable(sys::AbstractSystem, sym)
 end
 
 function SymbolicIndexingInterface.is_variable(sys::AbstractSystem, sym::Symbol)
+    if has_index_cache(sys) && (ic = get_index_cache(sys)) !== nothing
+        return haskey(ic.unknown_idx, hash(sym))
+    end
     return any(isequal(sym), getname.(variable_symbols(sys))) ||
            count('₊', string(sym)) == 1 &&
            count(isequal(sym), Symbol.(nameof(sys), :₊, getname.(variable_symbols(sys)))) ==
@@ -233,6 +398,9 @@ function SymbolicIndexingInterface.variable_index(sys::AbstractSystem, sym)
 end
 
 function SymbolicIndexingInterface.variable_index(sys::AbstractSystem, sym::Symbol)
+    if has_index_cache(sys) && (ic = get_index_cache(sys)) !== nothing
+        return get(ic.unknown_idx, h, nothing)
+    end
     idx = findfirst(isequal(sym), getname.(variable_symbols(sys)))
     if idx !== nothing
         return idx
@@ -257,12 +425,14 @@ function SymbolicIndexingInterface.is_parameter(sys::AbstractSystem, sym)
         ic = get_index_cache(sys)
         h = getsymbolhash(sym)
         return if haskey(ic.param_idx, h) || haskey(ic.discrete_idx, h) ||
-                  haskey(ic.constant_idx, h) || haskey(ic.dependent_idx, h)
+                  haskey(ic.constant_idx, h) || haskey(ic.dependent_idx, h) ||
+                  haskey(ic.nonnumeric_idx, h)
             true
         else
             h = getsymbolhash(default_toterm(sym))
             haskey(ic.param_idx, h) || haskey(ic.discrete_idx, h) ||
-                haskey(ic.constant_idx, h) || haskey(ic.dependent_idx, h)
+                haskey(ic.constant_idx, h) || haskey(ic.dependent_idx, h) ||
+                haskey(ic.nonnumeric_idx, h)
         end
     end
     return any(isequal(sym), parameter_symbols(sys)) ||
@@ -270,6 +440,9 @@ function SymbolicIndexingInterface.is_parameter(sys::AbstractSystem, sym)
 end
 
 function SymbolicIndexingInterface.is_parameter(sys::AbstractSystem, sym::Symbol)
+    if has_index_cache(sys) && (ic = get_index_cache(sys)) !== nothing
+        return ParameterIndex(ic, sym) !== nothing
+    end
     return any(isequal(sym), getname.(parameter_symbols(sys))) ||
            count('₊', string(sym)) == 1 &&
            count(isequal(sym),
@@ -282,7 +455,6 @@ function SymbolicIndexingInterface.parameter_index(sys::AbstractSystem, sym)
     end
     if has_index_cache(sys) && get_index_cache(sys) !== nothing
         ic = get_index_cache(sys)
-        h = getsymbolhash(sym)
         return if (idx = ParameterIndex(ic, sym)) !== nothing
             idx
         elseif (idx = ParameterIndex(ic, default_toterm(sym))) !== nothing
@@ -300,6 +472,9 @@ function SymbolicIndexingInterface.parameter_index(sys::AbstractSystem, sym)
 end
 
 function SymbolicIndexingInterface.parameter_index(sys::AbstractSystem, sym::Symbol)
+    if has_index_cache(sys) && (ic = get_index_cache(sys)) !== nothing
+        return ParameterIndex(ic, sym)
+    end
     idx = findfirst(isequal(sym), getname.(parameter_symbols(sys)))
     if idx !== nothing
         return idx
@@ -2113,6 +2288,10 @@ end
 
 keytype(::Type{<:Pair{T, V}}) where {T, V} = T
 function Symbolics.substitute(sys::AbstractSystem, rules::Union{Vector{<:Pair}, Dict})
+    if has_continuous_domain(sys) && get_continuous_events(sys) !== nothing ||
+       has_discrete_events(sys) && get_discrete_events(sys) !== nothing
+        @warn "`substitute` only supports performing substitutions in equations. This system has events, which will not be updated."
+    end
     if keytype(eltype(rules)) <: Symbol
         dict = todict(rules)
         systems = get_systems(sys)
@@ -2142,3 +2321,61 @@ function process_parameter_dependencies(pdeps, ps)
     end
     return pdeps, ps
 end
+
+"""
+    dump_parameters(sys::AbstractSystem)
+
+Return an array of `NamedTuple`s containing the metadata associated with each parameter in
+`sys`. Also includes the default value of the parameter, if provided.
+
+```@example
+using ModelingToolkit
+using DynamicQuantities
+using ModelingToolkit: t, D
+
+@parameters p = 1.0, [description = "My parameter", tunable = false] q = 2.0, [description = "Other parameter"]
+@variables x(t) = 3.0 [unit = u"m"]
+@named sys = ODESystem(Equation[], t, [x], [p, q])
+ModelingToolkit.dump_parameters(sys)
+```
+
+See also: [`ModelingToolkit.dump_variable_metadata`](@ref), [`ModelingToolkit.dump_unknowns`](@ref)
+"""
+function dump_parameters(sys::AbstractSystem)
+    defs = defaults(sys)
+    map(dump_variable_metadata.(parameters(sys))) do meta
+        if haskey(defs, meta.var)
+            meta = merge(meta, (; default = defs[meta.var]))
+        end
+        meta
+    end
+end
+
+"""
+    dump_unknowns(sys::AbstractSystem)
+
+Return an array of `NamedTuple`s containing the metadata associated with each unknown in
+`sys`. Also includes the default value of the unknown, if provided.
+
+```@example
+using ModelingToolkit
+using DynamicQuantities
+using ModelingToolkit: t, D
+
+@parameters p = 1.0, [description = "My parameter", tunable = false] q = 2.0, [description = "Other parameter"]
+@variables x(t) = 3.0 [unit = u"m"]
+@named sys = ODESystem(Equation[], t, [x], [p, q])
+ModelingToolkit.dump_unknowns(sys)
+```
+
+See also: [`ModelingToolkit.dump_variable_metadata`](@ref), [`ModelingToolkit.dump_parameters`](@ref)
+"""
+function dump_unknowns(sys::AbstractSystem)
+    defs = defaults(sys)
+    map(dump_variable_metadata.(unknowns(sys))) do meta
+        if haskey(defs, meta.var)
+            meta = merge(meta, (; default = defs[meta.var]))
+        end
+        meta
+    end
+end
diff --git a/src/systems/diffeqs/abstractodesystem.jl b/src/systems/diffeqs/abstractodesystem.jl
index 34899e3361..b6e8c50563 100644
--- a/src/systems/diffeqs/abstractodesystem.jl
+++ b/src/systems/diffeqs/abstractodesystem.jl
@@ -152,22 +152,8 @@ function generate_function(sys::AbstractODESystem, dvs = unknowns(sys),
         ddvs = implicit_dae ? map(Differential(get_iv(sys)), dvs) :
                nothing,
         isdde = false,
+        wrap_code = nothing,
         kwargs...)
-    array_vars = Dict{Any, Vector{Int}}()
-    for (j, x) in enumerate(dvs)
-        if istree(x) && operation(x) == getindex
-            arg = arguments(x)[1]
-            inds = get!(() -> Int[], array_vars, arg)
-            push!(inds, j)
-        end
-    end
-    subs = Dict()
-    for (k, inds) in array_vars
-        if inds == (inds′ = inds[1]:inds[end])
-            inds = inds′
-        end
-        subs[k] = term(view, Sym{Any}(Symbol("ˍ₋arg1")), inds)
-    end
     if isdde
         eqs = delay_to_function(sys)
     else
@@ -180,13 +166,15 @@ function generate_function(sys::AbstractODESystem, dvs = unknowns(sys),
     # substitute x(t) by just x
     rhss = implicit_dae ? [_iszero(eq.lhs) ? eq.rhs : eq.rhs - eq.lhs for eq in eqs] :
            [eq.rhs for eq in eqs]
-    rhss = fast_substitute(rhss, subs)
 
     # TODO: add an optional check on the ordering of observed equations
     u = map(x -> time_varying_as_func(value(x), sys), dvs)
     p = map.(x -> time_varying_as_func(value(x), sys), reorder_parameters(sys, ps))
     t = get_iv(sys)
 
+    if wrap_code === nothing
+        wrap_code = (identity, identity)
+    end
     if isdde
         build_function(rhss, u, DDE_HISTORY_FUN, p..., t; kwargs...)
     else
@@ -195,10 +183,12 @@ function generate_function(sys::AbstractODESystem, dvs = unknowns(sys),
         if implicit_dae
             build_function(rhss, ddvs, u, p..., t; postprocess_fbody = pre,
                 states = sol_states,
+                wrap_code = wrap_code .∘ wrap_array_vars(sys, rhss; dvs),
                 kwargs...)
         else
             build_function(rhss, u, p..., t; postprocess_fbody = pre,
                 states = sol_states,
+                wrap_code = wrap_code .∘ wrap_array_vars(sys, rhss; dvs),
                 kwargs...)
         end
     end
@@ -775,6 +765,16 @@ function get_u0_p(sys,
     if parammap !== nothing
         defs = mergedefaults(defs, parammap, ps)
     end
+    if u0map isa Vector && eltype(u0map) <: Pair
+        u0map = Dict(u0map)
+    end
+    if u0map isa Dict
+        allobs = Set(getproperty.(observed(sys), :lhs))
+        if any(in(allobs), keys(u0map))
+            u0s_in_obs = filter(in(allobs), keys(u0map))
+            @warn "Observed variables cannot assigned initial values. Initial values for $u0s_in_obs will be ignored."
+        end
+    end
     defs = mergedefaults(defs, u0map, dvs)
     for (k, v) in defs
         if Symbolics.isarraysymbolic(k)
diff --git a/src/systems/index_cache.jl b/src/systems/index_cache.jl
index 9228320c36..a7e6c8b157 100644
--- a/src/systems/index_cache.jl
+++ b/src/systems/index_cache.jl
@@ -40,10 +40,16 @@ function IndexCache(sys::AbstractSystem)
     let idx = 1
         for sym in unks
             h = getsymbolhash(sym)
-            if Symbolics.isarraysymbolic(sym)
-                unk_idxs[h] = idx:(idx + length(sym) - 1)
+            sym_idx = if Symbolics.isarraysymbolic(sym)
+                idx:(idx + length(sym) - 1)
             else
-                unk_idxs[h] = idx
+                idx
+            end
+            unk_idxs[h] = sym_idx
+
+            if hasname(sym)
+                h = hash(getname(sym))
+                unk_idxs[h] = sym_idx
             end
             idx += length(sym)
         end
@@ -122,6 +128,10 @@ function IndexCache(sys::AbstractSystem)
                 idxs[h] = (i, j)
                 h = getsymbolhash(default_toterm(p))
                 idxs[h] = (i, j)
+                if hasname(p)
+                    h = hash(getname(p))
+                    idxs[h] = (i, j)
+                end
             end
             push!(buffer_sizes, BufferTemplate(T, length(buf)))
         end
@@ -151,7 +161,7 @@ end
 
 function ParameterIndex(ic::IndexCache, p, sub_idx = ())
     p = unwrap(p)
-    h = getsymbolhash(p)
+    h = p isa Symbol ? hash(p) : getsymbolhash(p)
     return if haskey(ic.param_idx, h)
         ParameterIndex(SciMLStructures.Tunable(), (ic.param_idx[h]..., sub_idx...))
     elseif haskey(ic.discrete_idx, h)
diff --git a/src/systems/model_parsing.jl b/src/systems/model_parsing.jl
index a6bd0ce8c9..0f5e82aad6 100644
--- a/src/systems/model_parsing.jl
+++ b/src/systems/model_parsing.jl
@@ -120,9 +120,7 @@ function parse_variable_def!(dict, mod, arg, varclass, kwargs;
         (:misc, VariableMisc),
         (:disturbance, VariableDisturbance),
         (:tunable, VariableTunable),
-        (:dist, VariableDistribution),
-        (:binary, VariableBinary),
-        (:integer, VariableInteger)]
+        (:dist, VariableDistribution)]
 
     arg isa LineNumberNode && return
     MLStyle.@match arg begin
diff --git a/src/systems/optimization/optimizationsystem.jl b/src/systems/optimization/optimizationsystem.jl
index 12b44074e6..71b37ea978 100644
--- a/src/systems/optimization/optimizationsystem.jl
+++ b/src/systems/optimization/optimizationsystem.jl
@@ -258,8 +258,9 @@ function DiffEqBase.OptimizationProblem{iip}(sys::OptimizationSystem, u0map,
     if isnothing(lb) && isnothing(ub) # use the symbolically specified bounds
         lb = first.(getbounds.(dvs))
         ub = last.(getbounds.(dvs))
-        lb[isbinaryvar.(dvs)] .= 0
-        ub[isbinaryvar.(dvs)] .= 1
+        isboolean = symtype.(unwrap.(dvs)) .<: Bool
+        lb[isboolean] .= 0
+        ub[isboolean] .= 1
     else # use the user supplied variable bounds
         xor(isnothing(lb), isnothing(ub)) &&
             throw(ArgumentError("Expected both `lb` and `ub` to be supplied"))
@@ -269,7 +270,7 @@ function DiffEqBase.OptimizationProblem{iip}(sys::OptimizationSystem, u0map,
             throw(ArgumentError("Expected both `ub` to be of the same length as the vector of optimization variables"))
     end
 
-    int = isintegervar.(dvs) .| isbinaryvar.(dvs)
+    int = symtype.(unwrap.(dvs)) .<: Integer
 
     defs = defaults(sys)
     defs = mergedefaults(defs, parammap, ps)
@@ -476,8 +477,9 @@ function OptimizationProblemExpr{iip}(sys::OptimizationSystem, u0map,
     if isnothing(lb) && isnothing(ub) # use the symbolically specified bounds
         lb = first.(getbounds.(dvs))
         ub = last.(getbounds.(dvs))
-        lb[isbinaryvar.(dvs)] .= 0
-        ub[isbinaryvar.(dvs)] .= 1
+        isboolean = symtype.(unwrap.(dvs)) .<: Bool
+        lb[isboolean] .= 0
+        ub[isboolean] .= 1
     else # use the user supplied variable bounds
         xor(isnothing(lb), isnothing(ub)) &&
             throw(ArgumentError("Expected both `lb` and `ub` to be supplied"))
@@ -487,7 +489,7 @@ function OptimizationProblemExpr{iip}(sys::OptimizationSystem, u0map,
             throw(ArgumentError("Expected `ub` to be of the same length as the vector of optimization variables"))
     end
 
-    int = isintegervar.(dvs) .| isbinaryvar.(dvs)
+    int = symtype.(unwrap.(dvs)) .<: Integer
 
     defs = defaults(sys)
     defs = mergedefaults(defs, parammap, ps)
diff --git a/src/systems/parameter_buffer.jl b/src/systems/parameter_buffer.jl
index 897225c47f..62e3c39f72 100644
--- a/src/systems/parameter_buffer.jl
+++ b/src/systems/parameter_buffer.jl
@@ -128,7 +128,7 @@ end
 function split_into_buffers(raw::AbstractArray, buf; recurse = true)
     idx = 1
     function _helper(buf_v; recurse = true)
-        if eltype(buf_v) isa AbstractArray && recurse
+        if eltype(buf_v) <: AbstractArray && recurse
             return _helper.(buf_v; recurse = false)
         else
             res = reshape(raw[idx:(idx + length(buf_v) - 1)], size(buf_v))
@@ -139,6 +139,19 @@ function split_into_buffers(raw::AbstractArray, buf; recurse = true)
     return Tuple(_helper(buf_v; recurse) for buf_v in buf)
 end
 
+function update_tuple_of_buffers(raw::AbstractArray, buf)
+    idx = 1
+    function _helper(buf_v)
+        if eltype(buf_v) <: AbstractArray
+            _helper.(buf_v)
+        else
+            copyto!(buf_v, view(raw, idx:(idx + length(buf_v) - 1)))
+            idx += length(buf_v)
+        end
+    end
+    _helper.(buf)
+end
+
 SciMLStructures.isscimlstructure(::MTKParameters) = true
 
 SciMLStructures.ismutablescimlstructure(::MTKParameters) = true
@@ -151,7 +164,7 @@ for (Portion, field) in [(SciMLStructures.Tunable, :tunable)
         repack = let as_vector = as_vector, p = p
             function (new_val)
                 if new_val !== as_vector
-                    p.$field = split_into_buffers(new_val, p.$field)
+                    update_tuple_of_buffers(new_val, p.$field)
                 end
                 if p.dependent_update_iip !== nothing
                     p.dependent_update_iip(ArrayPartition(p.dependent), p...)
diff --git a/src/systems/unit_check.jl b/src/systems/unit_check.jl
index 6ff59f08fa..677d29895f 100644
--- a/src/systems/unit_check.jl
+++ b/src/systems/unit_check.jl
@@ -1,5 +1,5 @@
 #For dispatching get_unit
-const Conditional = Union{typeof(ifelse), typeof(IfElse.ifelse)}
+const Conditional = Union{typeof(ifelse)}
 const Comparison = Union{typeof.([==, !=, ≠, <, <=, ≤, >, >=, ≥])...}
 
 struct ValidationError <: Exception
diff --git a/src/systems/validation.jl b/src/systems/validation.jl
index f0c5fa95fd..ea458d8b7e 100644
--- a/src/systems/validation.jl
+++ b/src/systems/validation.jl
@@ -1,6 +1,6 @@
 module UnitfulUnitCheck
 
-using ..ModelingToolkit, Symbolics, SciMLBase, Unitful, IfElse, RecursiveArrayTools
+using ..ModelingToolkit, Symbolics, SciMLBase, Unitful, RecursiveArrayTools
 using ..ModelingToolkit: ValidationError,
                          ModelingToolkit, Connection, instream, JumpType, VariableUnit,
                          get_systems,
diff --git a/src/variables.jl b/src/variables.jl
index ca6d1b6954..9a9dbe4364 100644
--- a/src/variables.jl
+++ b/src/variables.jl
@@ -15,6 +15,69 @@ Symbolics.option_to_metadata_type(::Val{:irreducible}) = VariableIrreducible
 Symbolics.option_to_metadata_type(::Val{:state_priority}) = VariableStatePriority
 Symbolics.option_to_metadata_type(::Val{:misc}) = VariableMisc
 
+"""
+    dump_variable_metadata(var)
+
+Return all the metadata associated with symbolic variable `var` as a `NamedTuple`.
+
+```@example
+using ModelingToolkit
+
+@parameters p::Int [description = "My description", bounds = (0.5, 1.5)]
+ModelingToolkit.dump_variable_metadata(p)
+```
+"""
+function dump_variable_metadata(var)
+    uvar = unwrap(var)
+    vartype, name = get(uvar.metadata, VariableSource, (:unknown, :unknown))
+    shape = Symbolics.shape(var)
+    if shape == ()
+        shape = nothing
+    end
+    unit = get(uvar.metadata, VariableUnit, nothing)
+    connect = get(uvar.metadata, VariableConnectType, nothing)
+    noise = get(uvar.metadata, VariableNoiseType, nothing)
+    input = isinput(uvar) || nothing
+    output = isoutput(uvar) || nothing
+    irreducible = get(uvar.metadata, VariableIrreducible, nothing)
+    state_priority = get(uvar.metadata, VariableStatePriority, nothing)
+    misc = get(uvar.metadata, VariableMisc, nothing)
+    bounds = hasbounds(uvar) ? getbounds(uvar) : nothing
+    desc = getdescription(var)
+    if desc == ""
+        desc = nothing
+    end
+    guess = getguess(uvar)
+    disturbance = isdisturbance(uvar) || nothing
+    tunable = istunable(uvar, isparameter(uvar))
+    dist = getdist(uvar)
+    type = symtype(uvar)
+
+    meta = (
+        var = var,
+        vartype,
+        name,
+        shape,
+        unit,
+        connect,
+        noise,
+        input,
+        output,
+        irreducible,
+        state_priority,
+        misc,
+        bounds,
+        desc,
+        guess,
+        disturbance,
+        tunable,
+        dist,
+        type
+    )
+
+    return NamedTuple(k => v for (k, v) in pairs(meta) if v !== nothing)
+end
+
 abstract type AbstractConnectType end
 struct Equality <: AbstractConnectType end # Equality connection
 struct Flow <: AbstractConnectType end     # sum to 0
@@ -243,7 +306,7 @@ Symbolics.option_to_metadata_type(::Val{:tunable}) = VariableTunable
 istunable(x::Num, args...) = istunable(Symbolics.unwrap(x), args...)
 
 """
-    istunable(x, default = false)
+    istunable(x, default = true)
 
 Determine whether symbolic variable `x` is marked as a tunable for an automatic tuning algorithm.
 
@@ -257,7 +320,7 @@ Create a tunable parameter by
 
 See also [`tunable_parameters`](@ref), [`getbounds`](@ref)
 """
-function istunable(x, default = false)
+function istunable(x, default = true)
     p = Symbolics.getparent(x, nothing)
     p === nothing || (x = p)
     Symbolics.getmetadata(x, VariableTunable, default)
@@ -302,7 +365,7 @@ end
 ## System interface
 
 """
-    tunable_parameters(sys, p = parameters(sys); default=false)
+    tunable_parameters(sys, p = parameters(sys); default=true)
 
 Get all parameters of `sys` that are marked as `tunable`.
 
@@ -316,7 +379,7 @@ Create a tunable parameter by
 
 See also [`getbounds`](@ref), [`istunable`](@ref)
 """
-function tunable_parameters(sys, p = parameters(sys); default = false)
+function tunable_parameters(sys, p = parameters(sys); default = true)
     filter(x -> istunable(x, default), p)
 end
 
@@ -376,40 +439,6 @@ function hasdescription(x)
     getdescription(x) != ""
 end
 
-## binary variables =================================================================
-struct VariableBinary end
-Symbolics.option_to_metadata_type(::Val{:binary}) = VariableBinary
-
-isbinaryvar(x::Num) = isbinaryvar(Symbolics.unwrap(x))
-
-"""
-    isbinaryvar(x)
-
-Determine if a variable is binary.
-"""
-function isbinaryvar(x)
-    p = Symbolics.getparent(x, nothing)
-    p === nothing || (x = p)
-    return Symbolics.getmetadata(x, VariableBinary, false)
-end
-
-## integer variables =================================================================
-struct VariableInteger end
-Symbolics.option_to_metadata_type(::Val{:integer}) = VariableInteger
-
-isintegervar(x::Num) = isintegervar(Symbolics.unwrap(x))
-
-"""
-    isintegervar(x)
-
-Determine if a variable is an integer.
-"""
-function isintegervar(x)
-    p = Symbolics.getparent(x, nothing)
-    p === nothing || (x = p)
-    return Symbolics.getmetadata(x, VariableInteger, false)
-end
-
 ## Brownian
 """
     tobrownian(s::Sym)
diff --git a/test/domain_connectors.jl b/test/domain_connectors.jl
index 572f440681..9a43d2938f 100644
--- a/test/domain_connectors.jl
+++ b/test/domain_connectors.jl
@@ -1,7 +1,6 @@
 using ModelingToolkit
 using ModelingToolkit: t_nounits as t, D_nounits as D
 using Test
-using IfElse: ifelse
 
 @connector function HydraulicPort(; p_int, name)
     pars = @parameters begin
diff --git a/test/downstream/linearize.jl b/test/downstream/linearize.jl
index 302fabed94..c1aa6f6724 100644
--- a/test/downstream/linearize.jl
+++ b/test/downstream/linearize.jl
@@ -163,7 +163,7 @@ end
 function saturation(; y_max, y_min = y_max > 0 ? -y_max : -Inf, name)
     @variables u(t)=0 y(t)=0
     @parameters y_max=y_max y_min=y_min
-    ie = ModelingToolkit.IfElse.ifelse
+    ie = ifelse
     eqs = [
     # The equation below is equivalent to y ~ clamp(u, y_min, y_max)
         y ~ ie(u > y_max, y_max, ie((y_min < u) & (u < y_max), u, y_min))
diff --git a/test/dq_units.jl b/test/dq_units.jl
index a4bf6cba9b..75bbd4c4a9 100644
--- a/test/dq_units.jl
+++ b/test/dq_units.jl
@@ -1,4 +1,4 @@
-using ModelingToolkit, OrdinaryDiffEq, JumpProcesses, IfElse, DynamicQuantities
+using ModelingToolkit, OrdinaryDiffEq, JumpProcesses, DynamicQuantities
 using Test
 MT = ModelingToolkit
 using ModelingToolkit: t, D
diff --git a/test/generate_custom_function.jl b/test/generate_custom_function.jl
new file mode 100644
index 0000000000..9b64b20c12
--- /dev/null
+++ b/test/generate_custom_function.jl
@@ -0,0 +1,52 @@
+using ModelingToolkit
+using ModelingToolkit: t_nounits as t, D_nounits as D
+
+@variables x(t) y(t)[1:3]
+@parameters p1=1.0 p2[1:3]=[1.0, 2.0, 3.0] p3::Int=1 p4::Bool=false
+
+sys = complete(ODESystem(Equation[], t, [x; y], [p1, p2, p3, p4]; name = :sys))
+u0 = [1.0, 2.0, 3.0, 4.0]
+p = ModelingToolkit.MTKParameters(sys, [])
+
+fn1 = generate_custom_function(sys, x + y[1] + p1 + p2[1] + p3 * t; expression = Val(false))
+@test fn1(u0, p, 0.0) == 5.0
+
+fn2 = generate_custom_function(
+    sys, x + y[1] + p1 + p2[1] + p3 * t, [x], [p1, p2, p3]; expression = Val(false))
+@test fn1(u0, p, 0.0) == 5.0
+
+fn3_oop, fn3_iip = generate_custom_function(
+    sys, [x + y[2], y[3] + p2[2], p1 + p3, 3t]; expression = Val(false))
+
+buffer = zeros(4)
+fn3_iip(buffer, u0, p, 1.0)
+@test buffer == [4.0, 6.0, 2.0, 3.0]
+@test fn3_oop(u0, p, 1.0) == [4.0, 6.0, 2.0, 3.0]
+
+fn4 = generate_custom_function(sys, ifelse(p4, p1, p2[2]); expression = Val(false))
+@test fn4(u0, p, 1.0) == 2.0
+fn5 = generate_custom_function(sys, ifelse(!p4, p1, p2[2]); expression = Val(false))
+@test fn5(u0, p, 1.0) == 1.0
+
+@variables x y[1:3]
+sys = complete(NonlinearSystem(Equation[], [x; y], [p1, p2, p3, p4]; name = :sys))
+
+fn1 = generate_custom_function(sys, x + y[1] + p1 + p2[1] + p3; expression = Val(false))
+@test fn1(u0, p) == 6.0
+
+fn2 = generate_custom_function(
+    sys, x + y[1] + p1 + p2[1] + p3, [x], [p1, p2, p3]; expression = Val(false))
+@test fn1(u0, p) == 6.0
+
+fn3_oop, fn3_iip = generate_custom_function(
+    sys, [x + y[2], y[3] + p2[2], p1 + p3]; expression = Val(false))
+
+buffer = zeros(3)
+fn3_iip(buffer, u0, p)
+@test buffer == [4.0, 6.0, 2.0]
+@test fn3_oop(u0, p, 1.0) == [4.0, 6.0, 2.0]
+
+fn4 = generate_custom_function(sys, ifelse(p4, p1, p2[2]); expression = Val(false))
+@test fn4(u0, p, 1.0) == 2.0
+fn5 = generate_custom_function(sys, ifelse(!p4, p1, p2[2]); expression = Val(false))
+@test fn5(u0, p, 1.0) == 1.0
diff --git a/test/model_parsing.jl b/test/model_parsing.jl
index 46f455497e..70b90b669f 100644
--- a/test/model_parsing.jl
+++ b/test/model_parsing.jl
@@ -273,14 +273,14 @@ end
 
 @testset "Metadata in variables" begin
     metadata = Dict(:description => "Variable to test metadata in the Model.structure",
-        :input => true, :bounds => (-1, 1), :connection_type => :Flow, :integer => true,
-        :binary => false, :tunable => false, :disturbance => true, :dist => Normal(1, 1))
+        :input => true, :bounds => (-1, 1), :connection_type => :Flow,
+        :tunable => false, :disturbance => true, :dist => Normal(1, 1))
 
     @connector MockMeta begin
         m(t),
         [description = "Variable to test metadata in the Model.structure",
-            input = true, bounds = (-1, 1), connect = Flow, integer = true,
-            binary = false, tunable = false, disturbance = true, dist = Normal(1, 1)]
+            input = true, bounds = (-1, 1), connect = Flow,
+            tunable = false, disturbance = true, dist = Normal(1, 1)]
     end
 
     for (k, v) in metadata
diff --git a/test/runtests.jl b/test/runtests.jl
index 610ed55354..96bdbbf06e 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -63,6 +63,7 @@ end
         @safetestset "FuncAffect Test" include("funcaffect.jl")
         @safetestset "Constants Test" include("constants.jl")
         @safetestset "Parameter Dependency Test" include("parameter_dependencies.jl")
+        @safetestset "Generate Custom Function Test" include("generate_custom_function.jl")
     end
 
     if GROUP == "All" || GROUP == "InterfaceII"
diff --git a/test/state_selection.jl b/test/state_selection.jl
index 34e4767460..c1e4328715 100644
--- a/test/state_selection.jl
+++ b/test/state_selection.jl
@@ -1,4 +1,4 @@
-using ModelingToolkit, OrdinaryDiffEq, IfElse, Test
+using ModelingToolkit, OrdinaryDiffEq, Test
 using ModelingToolkit: t_nounits as t, D_nounits as D
 
 sts = @variables x1(t) x2(t) x3(t) x4(t)
@@ -249,7 +249,7 @@ let
     Δp2 = p2
 
     eqs = [+flow(xm, Δp1) ~ rho1 * dV1 + drho1 * V1
-           0 ~ IfElse.ifelse(w > 0.5,
+           0 ~ ifelse(w > 0.5,
                (0) - (rho2 * dV2 + drho2 * V2),
                (-flow(xm, Δp2)) - (rho2 * dV2 + drho2 * V2))
            V1 ~ (l_1f + w) * A_1f
@@ -265,7 +265,7 @@ let
            D(w) ~ dw
            D(dw) ~ ddw
            xf ~ 20e-3 * (1 - cos(2 * π * 5 * t))
-           0 ~ IfElse.ifelse(w > 0.5,
+           0 ~ ifelse(w > 0.5,
                (m_total * ddw) - (p1 * A_1f - p2 * A_2f - damp * dw),
                (m_total * ddw) - (p1 * A_1f - p2 * A_2f))]
     # ----------------------------------------------------------------------------
diff --git a/test/test_variable_metadata.jl b/test/test_variable_metadata.jl
index 886458302e..cdaeae8b93 100644
--- a/test/test_variable_metadata.jl
+++ b/test/test_variable_metadata.jl
@@ -4,31 +4,43 @@ using ModelingToolkit
 @variables u [bounds = (-1, 1)]
 @test getbounds(u) == (-1, 1)
 @test hasbounds(u)
+@test ModelingToolkit.dump_variable_metadata(u).bounds == (-1, 1)
 
 @variables y
 @test !hasbounds(y)
+@test !haskey(ModelingToolkit.dump_variable_metadata(y), :bounds)
 
 # Guess
 @variables y [guess = 0]
 @test getguess(y) === 0
 @test hasguess(y) === true
+@test ModelingToolkit.dump_variable_metadata(y).guess == 0
 
 @variables y
 @test hasguess(y) === false
+@test !haskey(ModelingToolkit.dump_variable_metadata(y), :guess)
 
 # Disturbance
 @variables u [disturbance = true]
 @test isdisturbance(u)
+@test ModelingToolkit.dump_variable_metadata(u).disturbance
 
 @variables y
 @test !isdisturbance(y)
+@test !haskey(ModelingToolkit.dump_variable_metadata(y), :disturbance)
 
 # Tunable
 @parameters u [tunable = true]
 @test istunable(u)
+@test ModelingToolkit.dump_variable_metadata(u).tunable
+
+@parameters u2 [tunable = false]
+@test !istunable(u2)
+@test !ModelingToolkit.dump_variable_metadata(u2).tunable
 
 @parameters y
-@test !istunable(y)
+@test istunable(y)
+@test ModelingToolkit.dump_variable_metadata(y).tunable
 
 # Distributions
 struct FakeNormal end
@@ -36,20 +48,28 @@ d = FakeNormal()
 @parameters u [dist = d]
 @test hasdist(u)
 @test getdist(u) == d
+@test ModelingToolkit.dump_variable_metadata(u).dist == d
 
 @parameters y
 @test !hasdist(y)
+@test !haskey(ModelingToolkit.dump_variable_metadata(y), :dist)
 
 ## System interface
 @parameters t
 Dₜ = Differential(t)
 @variables x(t)=0 [bounds = (-10, 10)] u(t)=0 [input = true] y(t)=0 [output = true]
-@parameters T [tunable = true, bounds = (0, Inf)]
+@parameters T [bounds = (0, Inf)]
 @parameters k [tunable = true, bounds = (0, Inf)]
-@parameters k2
+@parameters k2 [tunable = false]
 eqs = [Dₜ(x) ~ (-k2 * x + k * u) / T
        y ~ x]
 sys = ODESystem(eqs, t, name = :tunable_first_order)
+unk_meta = ModelingToolkit.dump_unknowns(sys)
+@test length(unk_meta) == 3
+@test all(iszero, meta.default for meta in unk_meta)
+param_meta = ModelingToolkit.dump_parameters(sys)
+@test length(param_meta) == 3
+@test all(!haskey(meta, :default) for meta in param_meta)
 
 p = tunable_parameters(sys)
 sp = Set(p)
@@ -68,21 +88,23 @@ b = getbounds(sys)
 b = getbounds(sys, unknowns(sys))
 @test b[x] == (-10, 10)
 
-p = tunable_parameters(sys, default = true)
+p = tunable_parameters(sys, default = false)
 sp = Set(p)
 @test k ∈ sp
-@test T ∈ sp
-@test k2 ∈ sp
-@test length(p) == 3
+@test T ∉ sp
+@test k2 ∉ sp
+@test length(p) == 1
 
 ## Descriptions
 @variables u [description = "This is my input"]
 @test getdescription(u) == "This is my input"
 @test hasdescription(u)
+@test ModelingToolkit.dump_variable_metadata(u).desc == "This is my input"
 
 @variables u
 @test getdescription(u) == ""
 @test !hasdescription(u)
+@test !haskey(ModelingToolkit.dump_variable_metadata(u), :desc)
 
 @parameters t
 @variables u(t) [description = "A short description of u"]
@@ -90,19 +112,3 @@ sp = Set(p)
 @named sys = ODESystem([u ~ p], t)
 
 @test_nowarn show(stdout, "text/plain", sys)
-
-@testset "binary" begin
-    @parameters t
-    @variables u(t) [binary = true]
-    @parameters p [binary = true]
-    @test isbinaryvar(u)
-    @test isbinaryvar(p)
-end
-
-@testset "integer" begin
-    @parameters t
-    @variables u(t) [integer = true]
-    @parameters p [integer = true]
-    @test isintegervar(u)
-    @test isintegervar(p)
-end
diff --git a/test/units.jl b/test/units.jl
index fd737cf816..e170540270 100644
--- a/test/units.jl
+++ b/test/units.jl
@@ -1,4 +1,4 @@
-using ModelingToolkit, OrdinaryDiffEq, JumpProcesses, IfElse, Unitful
+using ModelingToolkit, OrdinaryDiffEq, JumpProcesses, Unitful
 using Test
 MT = ModelingToolkit
 UMT = ModelingToolkit.UnitfulUnitCheck
@@ -33,7 +33,7 @@ D = Differential(t)
 @test UMT.get_unit(t / τ) == UMT.unitless
 @test UMT.equivalent(UMT.get_unit(P - E / τ), u"MW")
 @test UMT.equivalent(UMT.get_unit(D(D(E))), u"MW/ms")
-@test UMT.get_unit(IfElse.ifelse(t > t, P, E / τ)) == u"MW"
+@test UMT.get_unit(ifelse(t > t, P, E / τ)) == u"MW"
 @test UMT.get_unit(1.0^(t / τ)) == UMT.unitless
 @test UMT.get_unit(exp(t / τ)) == UMT.unitless
 @test UMT.get_unit(sin(t / τ)) == UMT.unitless