diff --git a/.github/workflows/Tests.yml b/.github/workflows/Tests.yml
index a95f19a085..52c5482970 100644
--- a/.github/workflows/Tests.yml
+++ b/.github/workflows/Tests.yml
@@ -38,6 +38,7 @@ jobs:
           - Extensions
           - Downstream
           - RegressionI
+          - FMI
     uses: "SciML/.github/.github/workflows/tests.yml@v1"
     with:
       julia-version: "${{ matrix.version }}"
diff --git a/Project.toml b/Project.toml
index 0a6238685a..e4b2238170 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 = "9.61.0"
+version = "9.62.0"
 
 [deps]
 AbstractTrees = "1520ce14-60c1-5f80-bbc7-55ef81b5835c"
@@ -64,6 +64,7 @@ Unitful = "1986cc42-f94f-5a68-af5c-568840ba703d"
 BifurcationKit = "0f109fa4-8a5d-4b75-95aa-f515264e7665"
 ChainRulesCore = "d360d2e6-b24c-11e9-a2a3-2a2ae2dbcce4"
 DeepDiffs = "ab62b9b5-e342-54a8-a765-a90f495de1a6"
+FMI = "14a09403-18e3-468f-ad8a-74f8dda2d9ac"
 HomotopyContinuation = "f213a82b-91d6-5c5d-acf7-10f1c761b327"
 InfiniteOpt = "20393b10-9daf-11e9-18c9-8db751c92c57"
 LabelledArrays = "2ee39098-c373-598a-b85f-a56591580800"
@@ -72,6 +73,7 @@ LabelledArrays = "2ee39098-c373-598a-b85f-a56591580800"
 MTKBifurcationKitExt = "BifurcationKit"
 MTKChainRulesCoreExt = "ChainRulesCore"
 MTKDeepDiffsExt = "DeepDiffs"
+MTKFMIExt = "FMI"
 MTKHomotopyContinuationExt = "HomotopyContinuation"
 MTKInfiniteOptExt = "InfiniteOpt"
 MTKLabelledArraysExt = "LabelledArrays"
@@ -106,6 +108,7 @@ FindFirstFunctions = "1"
 ForwardDiff = "0.10.3"
 FunctionWrappers = "1.1"
 FunctionWrappersWrappers = "0.1"
+FMI = "0.14"
 Graphs = "1.5.2"
 HomotopyContinuation = "2.11"
 InfiniteOpt = "0.5"
@@ -132,7 +135,7 @@ RecursiveArrayTools = "3.26"
 Reexport = "0.2, 1"
 RuntimeGeneratedFunctions = "0.5.9"
 SCCNonlinearSolve = "1.0.0"
-SciMLBase = "2.71.1"
+SciMLBase = "2.72.1"
 SciMLStructures = "1.0"
 Serialization = "1"
 Setfield = "0.7, 0.8, 1"
@@ -142,9 +145,9 @@ SpecialFunctions = "0.7, 0.8, 0.9, 0.10, 1.0, 2"
 StaticArrays = "0.10, 0.11, 0.12, 1.0"
 StochasticDiffEq = "6.72.1"
 StochasticDelayDiffEq = "1.8.1"
-SymbolicIndexingInterface = "0.3.36"
-SymbolicUtils = "3.10"
-Symbolics = "6.22.1"
+SymbolicIndexingInterface = "0.3.37"
+SymbolicUtils = "3.10.1"
+Symbolics = "6.25"
 URIs = "1"
 UnPack = "0.1, 1.0"
 Unitful = "1.1"
diff --git a/ext/MTKFMIExt.jl b/ext/MTKFMIExt.jl
new file mode 100644
index 0000000000..5cfe9a82ef
--- /dev/null
+++ b/ext/MTKFMIExt.jl
@@ -0,0 +1,933 @@
+module MTKFMIExt
+
+using ModelingToolkit
+using SymbolicIndexingInterface
+using ModelingToolkit: t_nounits as t, D_nounits as D
+using DocStringExtensions
+import ModelingToolkit as MTK
+import SciMLBase
+import FMI
+
+"""
+    $(TYPEDSIGNATURES)
+
+A utility macro for FMI.jl functions that return a status. Will terminate on
+fatal statuses. Must be used as `@statuscheck FMI.fmiXFunction(...)` where
+`X` should be `2` or `3`. Has an edge case for handling tuples for
+`FMI.fmi2CompletedIntegratorStep`.
+"""
+macro statuscheck(expr)
+    @assert Meta.isexpr(expr, :call)
+    fn = expr.args[1]
+    @assert Meta.isexpr(fn, :.)
+    @assert fn.args[1] == :FMI
+    fnname = fn.args[2]
+
+    instance = expr.args[2]
+    is_v2 = startswith("fmi2", string(fnname))
+
+    fmiTrue = is_v2 ? FMI.fmi2True : FMI.fmi3True
+    fmiStatusOK = is_v2 ? FMI.fmi2StatusOK : FMI.fmi3StatusOK
+    fmiStatusWarning = is_v2 ? FMI.fmi2StatusWarning : FMI.fmi3StatusWarning
+    fmiStatusFatal = is_v2 ? FMI.fmi2StatusFatal : FMI.fmi3StatusFatal
+    fmiTerminate = is_v2 ? FMI.fmi2Terminate : FMI.fmi3Terminate
+    fmiFreeInstance! = is_v2 ? FMI.fmi2FreeInstance! : FMI.fmi3FreeInstance!
+    return quote
+        status = $expr
+        fnname = $fnname
+        if status !== nothing && ((status isa Tuple && status[1] == $fmiTrue) ||
+            (!(status isa Tuple) && status != $fmiStatusOK &&
+             status != $fmiStatusWarning))
+            if status != $fmiStatusFatal
+                $fmiTerminate(wrapper.instance)
+            end
+            $fmiFreeInstance!(wrapper.instance)
+            wrapper.instance = nothing
+            error("FMU Error in $fnname: status $status")
+        end
+    end |> esc
+end
+
+@static if !hasmethod(FMI.getValueReferencesAndNames, Tuple{FMI.fmi3ModelDescription})
+    """
+        $(TYPEDSIGNATURES)
+
+        This is type piracy, but FMI.jl is missing this implementation. It allows
+        `FMI.getStateValueReferencesAndNames` to work.
+    """
+    function FMI.getValueReferencesAndNames(
+            md::FMI.fmi3ModelDescription; vrs = md.valueReferences)
+        dict = Dict{FMI.fmi3ValueReference, Array{String}}()
+        for vr in vrs
+            dict[vr] = FMI.valueReferenceToString(md, vr)
+        end
+        return dict
+    end
+end
+
+"""
+    $(TYPEDSIGNATURES)
+
+A component that wraps an FMU loaded via FMI.jl. The FMI version (2 or 3) should be
+provided as a `Val` to the function. Supports Model Exchange and CoSimulation FMUs.
+All inputs, continuous variables and outputs must be `FMI.fmi2Real` or `FMI.fmi3Float64`.
+Does not support events or discrete variables in the FMU. Does not support automatic
+differentiation. Parameters of the FMU will have defaults corresponding to their initial
+values in the FMU specification. All other variables will not have a default. Hierarchical
+names in the FMU of the form `namespace.variable` are transformed into symbolic variables
+with the name `namespace__variable`.
+
+# Keyword Arguments
+
+- `fmu`: The FMU loaded via `FMI.loadFMU`.
+- `tolerance`: The tolerance to provide to the FMU. Not used for v3 FMUs since it is not
+  supported by FMI.jl.
+- `communication_step_size`: The periodic interval at which communication with CoSimulation
+  FMUs will occur. Must be provided for CoSimulation FMU components.
+- `reinitializealg`: The DAE initialization algorithm to use for the callback managing the
+  FMU. For CoSimulation FMUs whose states/outputs are used in algebraic equations of the
+  system, this needs to be an algorithm that will solve for the new algebraic variables.
+  For example, `OrdinaryDiffEqCore.BrownFullBasicInit()`.
+- `type`: Either `:ME` or `:CS` depending on whether `fmu` is a Model Exchange or
+  CoSimulation FMU respectively.
+- `name`: The name of the system.
+"""
+function MTK.FMIComponent(::Val{Ver}; fmu = nothing, tolerance = 1e-6,
+        communication_step_size = nothing, reinitializealg = SciMLBase.NoInit(), type, name) where {Ver}
+    if Ver != 2 && Ver != 3
+        throw(ArgumentError("FMI Version must be `2` or `3`"))
+    end
+    if type == :CS && communication_step_size === nothing
+        throw(ArgumentError("`communication_step_size` must be specified for Co-Simulation FMUs."))
+    end
+    # mapping from MTK variable to value reference
+    value_references = Dict()
+    # defaults
+    defs = Dict()
+    # unknowns of the system
+    states = []
+    # differential variables of the system
+    # this is a subset of `states` in the case where the FMU has multiple names for
+    # the same value reference.
+    diffvars = []
+    # variables that are derivatives of diffvars
+    dervars = []
+    # observed equations
+    observed = Equation[]
+    # need to separate observed equations for duplicate derivative variables
+    # since they aren't included in CS FMUs
+    der_observed = Equation[]
+
+    # parse states
+    fmi_variables_to_mtk_variables!(fmu, FMI.getStateValueReferencesAndNames(fmu),
+        value_references, diffvars, states, observed)
+    # create a symbolic variable __mtk_internal_u to pass to the relevant registered
+    # functions as the state vector
+    if isempty(diffvars)
+        # no differential variables
+        __mtk_internal_u = Float64[]
+    elseif type == :ME
+        # to avoid running into `structural_simplify` warnings about array variables
+        # and some unfortunate circular dependency issues, ME FMUs use an array of
+        # symbolics instead. This is also not worse off in performance
+        # because the former approach would allocate anyway.
+        # TODO: Can we avoid an allocation here using static arrays?
+        __mtk_internal_u = copy(diffvars)
+    elseif type == :CS
+        # CS FMUs do their own independent integration in a periodic callback, so their
+        # unknowns are discrete variables in the `ODESystem`. A default of `missing` allows
+        # them to be solved for during initialization.
+        @parameters __mtk_internal_u(t)[1:length(diffvars)]=missing [guess = diffvars]
+        push!(observed, __mtk_internal_u ~ copy(diffvars))
+    end
+
+    # parse derivatives of states
+    # the variables passed to `postprocess_variable` haven't been differentiated yet, so they
+    # should match one variable in states. That's the one this is the derivative of, and we
+    # keep track of this ordering
+    derivative_order = []
+    function derivative_order_postprocess(var)
+        idx = findfirst(isequal(var), states)
+        idx === nothing || push!(derivative_order, states[idx])
+        return var
+    end
+    fmi_variables_to_mtk_variables!(
+        fmu, FMI.getDerivateValueReferencesAndNames(fmu), value_references, dervars,
+        states, der_observed; postprocess_variable = derivative_order_postprocess)
+    @assert length(derivative_order) == length(dervars)
+
+    # parse the inputs to the FMU
+    inputs = []
+    fmi_variables_to_mtk_variables!(fmu, FMI.getInputValueReferencesAndNames(fmu),
+        value_references, inputs, states, observed; postprocess_variable = v -> MTK.setinput(
+            v, true))
+    # create a symbolic variable for the input buffer
+    __mtk_internal_x = copy(inputs)
+    if isempty(__mtk_internal_x)
+        __mtk_internal_x = Float64[]
+    end
+
+    # parse the outputs of the FMU
+    outputs = []
+    fmi_variables_to_mtk_variables!(fmu, FMI.getOutputValueReferencesAndNames(fmu),
+        value_references, outputs, states, observed; postprocess_variable = v -> MTK.setoutput(
+            v, true))
+    # create the output buffer. This is only required for CoSimulation to pass it to
+    # the callback affect
+    if type == :CS
+        if isempty(outputs)
+            __mtk_internal_o = Float64[]
+        else
+            @parameters __mtk_internal_o(t)[1:length(outputs)]=missing [guess = zeros(length(outputs))]
+            push!(observed, __mtk_internal_o ~ outputs)
+        end
+    end
+
+    # parse the parameters
+    params = []
+    # multiple names for the same parameter are treated as parameter dependencies.
+    parameter_dependencies = Equation[]
+    fmi_variables_to_mtk_variables!(
+        fmu, FMI.getParameterValueReferencesAndNames(fmu), value_references,
+        params, [], parameter_dependencies, defs; parameters = true)
+    # create a symbolic variable for the parameter buffer
+    __mtk_internal_p = copy(params)
+    if isempty(__mtk_internal_p)
+        __mtk_internal_p = Float64[]
+    end
+
+    derivative_value_references = UInt32[value_references[var] for var in dervars]
+    state_value_references = UInt32[value_references[var] for var in diffvars]
+    output_value_references = UInt32[value_references[var] for var in outputs]
+    input_value_references = UInt32[value_references[var] for var in inputs]
+    param_value_references = UInt32[value_references[var] for var in params]
+
+    # create a parameter for the instance wrapper
+    # this manages the creation and deallocation of FMU instances
+    buffer_length = length(diffvars) + length(outputs)
+    if Ver == 2
+        @parameters (wrapper::FMI2InstanceWrapper)(..)[1:buffer_length] = FMI2InstanceWrapper(
+            fmu, derivative_value_references, state_value_references, output_value_references,
+            param_value_references, input_value_references, tolerance)
+    else
+        @parameters (wrapper::FMI3InstanceWrapper)(..)[1:buffer_length] = FMI3InstanceWrapper(
+            fmu, derivative_value_references, state_value_references,
+            output_value_references, param_value_references, input_value_references)
+    end
+
+    # any additional initialization equations for the system
+    initialization_eqs = Equation[]
+
+    if type == :ME
+        # the wrapper is a callable struct which returns the state derivative and
+        # output values
+        # symbolic expression for calling the wrapper
+        call_expr = wrapper(__mtk_internal_u, __mtk_internal_x, __mtk_internal_p, t)
+
+        # differential and observed equations
+        diffeqs = Equation[]
+        for (i, var) in enumerate([dervars; outputs])
+            push!(diffeqs, var ~ call_expr[i])
+        end
+        for (var, dervar) in zip(derivative_order, dervars)
+            push!(diffeqs, D(var) ~ dervar)
+        end
+
+        # instance management callback which deallocates the instance when
+        # necessary and notifies the FMU of completed integrator steps
+        finalize_affect = MTK.FunctionalAffect(fmiFinalize!, [], [wrapper], [])
+        step_affect = MTK.FunctionalAffect(Returns(nothing), [], [], [])
+        instance_management_callback = MTK.SymbolicDiscreteCallback(
+            (t != t - 1), step_affect; finalize = finalize_affect, reinitializealg = reinitializealg)
+
+        push!(params, wrapper)
+        append!(observed, der_observed)
+    elseif type == :CS
+        _functor = if Ver == 2
+            FMI2CSFunctor(state_value_references, output_value_references)
+        else
+            FMI3CSFunctor(state_value_references, output_value_references)
+        end
+        @parameters (functor::(typeof(_functor)))(..)[1:(length(__mtk_internal_u) + length(__mtk_internal_o))] = _functor
+        # for co-simulation, we need to ensure the output buffer is solved for
+        # during initialization
+        for (i, x) in enumerate(collect(__mtk_internal_o))
+            push!(initialization_eqs,
+                x ~ functor(
+                    wrapper, __mtk_internal_u, __mtk_internal_x, __mtk_internal_p, t)[i])
+        end
+
+        diffeqs = Equation[]
+
+        # use `ImperativeAffect` for instance management here
+        cb_observed = (; inputs = __mtk_internal_x, params = copy(params),
+            t, wrapper, dt = communication_step_size)
+        cb_modified = (;)
+        # modify the outputs if present
+        if symbolic_type(__mtk_internal_o) != NotSymbolic()
+            cb_modified = (cb_modified..., outputs = __mtk_internal_o)
+        end
+        # modify the continuous state if present
+        if symbolic_type(__mtk_internal_u) != NotSymbolic()
+            cb_modified = (cb_modified..., states = __mtk_internal_u)
+        end
+        initialize_affect = MTK.ImperativeAffect(fmiCSInitialize!; observed = cb_observed,
+            modified = cb_modified, ctx = _functor)
+        finalize_affect = MTK.FunctionalAffect(fmiFinalize!, [], [wrapper], [])
+        # the callback affect performs the stepping
+        step_affect = MTK.ImperativeAffect(
+            fmiCSStep!; observed = cb_observed, modified = cb_modified, ctx = _functor)
+        instance_management_callback = MTK.SymbolicDiscreteCallback(
+            communication_step_size, step_affect; initialize = initialize_affect,
+            finalize = finalize_affect, reinitializealg = reinitializealg
+        )
+
+        # guarded in case there are no outputs/states and the variable is `[]`.
+        symbolic_type(__mtk_internal_o) == NotSymbolic() || push!(params, __mtk_internal_o)
+        symbolic_type(__mtk_internal_u) == NotSymbolic() || push!(params, __mtk_internal_u)
+
+        push!(params, wrapper, functor)
+    end
+
+    eqs = [observed; diffeqs]
+    return ODESystem(eqs, t, states, params; parameter_dependencies, defaults = defs,
+        discrete_events = [instance_management_callback], name, initialization_eqs)
+end
+
+"""
+    $(TYPEDSIGNATURES)
+
+A utility function which accepts an FMU `fmu` and a mapping from value reference to a
+list of associated names `varmap`. A symbolic variable is created for each name. The
+associated value reference is kept track of in `value_references`. In case there are
+multiple names for a value reference, the symbolic variable for the first name is pushed
+to `truevars`. All of the created symbolic variables are pushed to `allvars`. Observed
+equations equating identical variables are pushed to `obseqs`. `defs` is a dictionary of
+defaults.
+
+# Keyword Arguments
+- `parameters`: A boolean indicating whether to use `@parameters` for the symbolic
+  variables instead of `@variables`.
+- `postprocess_variable`: A function applied to each created variable that should
+  return the updated variable. This is useful to add metadata to variables.
+"""
+function fmi_variables_to_mtk_variables!(
+        fmu::Union{FMI.FMU2, FMI.FMU3}, varmap::AbstractDict,
+        value_references::AbstractDict, truevars, allvars,
+        obseqs, defs = Dict(); parameters = false, postprocess_variable = identity)
+    for (valRef, varnames) in varmap
+        stateT = FMI.dataTypeForValueReference(fmu, valRef)
+        snames = Symbol[]
+        ders = Int[]
+        for name in varnames
+            sname, der = parseFMIVariableName(name)
+            push!(snames, sname)
+            push!(ders, der)
+        end
+        if parameters
+            vars = [postprocess_variable(MTK.unwrap(only(@parameters $sname::stateT)))
+                    for sname in snames]
+        else
+            vars = [postprocess_variable(MTK.unwrap(only(@variables $sname(t)::stateT)))
+                    for sname in snames]
+        end
+        for i in eachindex(vars)
+            der = ders[i]
+            vars[i] = MTK.unwrap(vars[i])
+            for j in 1:der
+                vars[i] = D(vars[i])
+            end
+            vars[i] = MTK.default_toterm(vars[i])
+        end
+        for i in eachindex(vars)
+            if i == 1
+                push!(truevars, vars[i])
+            else
+                push!(obseqs, vars[i] ~ vars[1])
+            end
+            value_references[vars[i]] = valRef
+        end
+        append!(allvars, vars)
+        defval = FMI.getStartValue(fmu, valRef)
+        defs[vars[1]] = defval
+    end
+end
+
+"""
+    $(TYPEDSIGNATURES)
+
+Parse the string name of an FMI variable into a `Symbol` name for the corresponding
+MTK variable. Return the `Symbol` name and the number of times it is differentiated.
+"""
+function parseFMIVariableName(name::AbstractString)
+    name = replace(name, "." => "__")
+    der = 0
+    if startswith(name, "der(")
+        idx = findfirst(',', name)
+        if idx === nothing
+            name = @view name[5:(end - 1)]
+            der = 1
+        else
+            der = parse(Int, @view name[(idx + 1):(end - 1)])
+            name = @view name[5:(idx - 1)]
+        end
+    end
+    return Symbol(name), der
+end
+
+"""
+    $(TYPEDEF)
+
+A struct which manages instance creation and deallocation for v2 FMUs.
+
+# Fields
+
+$(TYPEDFIELDS)
+"""
+mutable struct FMI2InstanceWrapper
+    """
+    The FMU from `FMI.loadFMU`.
+    """
+    const fmu::FMI.FMU2
+    """
+    The value references for derivatives of states of the FMU, in the order that the
+    caller expects them to be returned when calling this struct.
+    """
+    const derivative_value_references::Vector{FMI.fmi2ValueReference}
+    """
+    The value references for the states of the FMU.
+    """
+    const state_value_references::Vector{FMI.fmi2ValueReference}
+    """
+    The value references for outputs of the FMU, in the order that the caller expects
+    them to be returned when calling this struct.
+    """
+    const output_value_references::Vector{FMI.fmi2ValueReference}
+    """
+    The parameter value references. These should be in the same order as the parameter
+    vector passed to functions involving this wrapper.
+    """
+    const param_value_references::Vector{FMI.fmi2ValueReference}
+    """
+    The input value references. These should be in the same order as the inputs passed
+    to functions involving this wrapper.
+    """
+    const input_value_references::Vector{FMI.fmi2ValueReference}
+    """
+    The tolerance with which to setup the FMU instance.
+    """
+    const tolerance::FMI.fmi2Real
+    """
+    The FMU instance, if present, and `nothing` otherwise.
+    """
+    instance::Union{FMI.FMU2Component{FMI.FMU2}, Nothing}
+end
+
+"""
+    $(TYPEDSIGNATURES)
+
+Create an `FMI2InstanceWrapper` with no instance.
+"""
+function FMI2InstanceWrapper(fmu, ders, states, outputs, params, inputs, tolerance)
+    FMI2InstanceWrapper(fmu, ders, states, outputs, params, inputs, tolerance, nothing)
+end
+
+Base.nameof(::FMI2InstanceWrapper) = :FMI2InstanceWrapper
+
+"""
+    $(TYPEDSIGNATURES)
+
+Common functionality for creating an instance of a v2 FMU. Does not check if
+`wrapper.instance` is already present, and overwrites the existing value with
+a new instance. `inputs` should be in the order of `wrapper.input_value_references`.
+`params` should be in the order of `wrapper.param_value_references`. `t` is the current
+time. Returns the created instance, which is also stored in `wrapper.instance`.
+"""
+function get_instance_common!(wrapper::FMI2InstanceWrapper, inputs, params, t)
+    wrapper.instance = FMI.fmi2Instantiate!(wrapper.fmu)::FMI.FMU2Component
+    if !isempty(inputs)
+        @statuscheck FMI.fmi2SetReal(wrapper.instance, wrapper.input_value_references,
+            Csize_t(length(wrapper.param_value_references)), inputs)
+    end
+    if !isempty(params)
+        @statuscheck FMI.fmi2SetReal(wrapper.instance, wrapper.param_value_references,
+            Csize_t(length(wrapper.param_value_references)), params)
+    end
+    @statuscheck FMI.fmi2SetupExperiment(
+        wrapper.instance, FMI.fmi2True, wrapper.tolerance, t, FMI.fmi2False, t)
+    @statuscheck FMI.fmi2EnterInitializationMode(wrapper.instance)
+    return wrapper.instance
+end
+
+"""
+    $(TYPEDSIGNATURES)
+
+Create an instance of a Model Exchange FMU. Use the existing instance in `wrapper` if
+present and create a new one otherwise. Return the instance.
+
+See `get_instance_common!` for a description of the arguments.
+"""
+function get_instance_ME!(wrapper::FMI2InstanceWrapper, inputs, params, t)
+    if wrapper.instance === nothing
+        get_instance_common!(wrapper, inputs, params, t)
+        @statuscheck FMI.fmi2ExitInitializationMode(wrapper.instance)
+        eventInfo = FMI.fmi2NewDiscreteStates(wrapper.instance)
+        @assert eventInfo.newDiscreteStatesNeeded == FMI.fmi2False
+        # TODO: Support FMU events
+        @statuscheck FMI.fmi2EnterContinuousTimeMode(wrapper.instance)
+    end
+
+    return wrapper.instance
+end
+
+"""
+    $(TYPEDSIGNATURES)
+
+Create an instance of a CoSimulation FMU. Use the existing instance in `wrapper` if
+present and create a new one otherwise. Return the instance.
+
+See `get_instance_common!` for a description of the arguments.
+"""
+function get_instance_CS!(wrapper::FMI2InstanceWrapper, states, inputs, params, t)
+    if wrapper.instance === nothing
+        get_instance_common!(wrapper, inputs, params, t)
+        if !isempty(states)
+            @statuscheck FMI.fmi2SetReal(wrapper.instance, wrapper.state_value_references,
+                Csize_t(length(wrapper.state_value_references)), states)
+        end
+        @statuscheck FMI.fmi2ExitInitializationMode(wrapper.instance)
+    end
+    return wrapper.instance
+end
+
+"""
+    $(TYPEDSIGNATURES)
+
+Call `fmiXCompletedIntegratorStep` with `noSetFMUStatePriorToCurrentPoint` as false.
+"""
+function partiallyCompleteIntegratorStep(wrapper::FMI2InstanceWrapper)
+    @statuscheck FMI.fmi2CompletedIntegratorStep(wrapper.instance, FMI.fmi2False)
+end
+
+"""
+    $(TYPEDSIGNATURES)
+
+If `wrapper.instance !== nothing`, terminate and free the instance. Also set
+`wrapper.instance` to `nothing`.
+"""
+function reset_instance!(wrapper::FMI2InstanceWrapper)
+    wrapper.instance === nothing && return
+    FMI.fmi2Terminate(wrapper.instance)
+    FMI.fmi2FreeInstance!(wrapper.instance)
+    wrapper.instance = nothing
+end
+
+"""
+    $(TYPEDEF)
+
+A struct which manages instance creation and deallocation for v3 FMUs.
+
+# Fields
+
+$(TYPEDFIELDS)
+"""
+mutable struct FMI3InstanceWrapper
+    """
+    The FMU from `FMI.loadFMU`.
+    """
+    const fmu::FMI.FMU3
+    """
+    The value references for derivatives of states of the FMU, in the order that the
+    caller expects them to be returned when calling this struct.
+    """
+    const derivative_value_references::Vector{FMI.fmi3ValueReference}
+    const state_value_references::Vector{FMI.fmi3ValueReference}
+    """
+    The value references for outputs of the FMU, in the order that the caller expects
+    them to be returned when calling this struct.
+    """
+    const output_value_references::Vector{FMI.fmi3ValueReference}
+    """
+    The parameter value references. These should be in the same order as the parameter
+    vector passed to functions involving this wrapper.
+    """
+    const param_value_references::Vector{FMI.fmi3ValueReference}
+    """
+    The input value references. These should be in the same order as the inputs passed
+    to functions involving this wrapper.
+    """
+    const input_value_references::Vector{FMI.fmi3ValueReference}
+    """
+    The FMU instance, if present, and `nothing` otherwise.
+    """
+    instance::Union{FMI.FMU3Instance{FMI.FMU3}, Nothing}
+end
+
+"""
+    $(TYPEDSIGNATURES)
+
+Create an `FMI3InstanceWrapper` with no instance.
+"""
+function FMI3InstanceWrapper(fmu, ders, states, outputs, params, inputs)
+    FMI3InstanceWrapper(fmu, ders, states, outputs, params, inputs, nothing)
+end
+
+Base.nameof(::FMI3InstanceWrapper) = :FMI3InstanceWrapper
+
+"""
+    $(TYPEDSIGNATURES)
+
+Common functionality for creating an instance of a v3 FMU. Since v3 FMUs need to be
+instantiated differently depending on the type, this assumes `wrapper.instance` is a
+freshly instantiated FMU which needs to be initialized. `inputs` should be in the order
+of `wrapper.input_value_references`. `params` should be in the order of
+`wrapper.param_value_references`. `t` is the current time. Returns `wrapper.instance`.
+"""
+function get_instance_common!(wrapper::FMI3InstanceWrapper, inputs, params, t)
+    if !isempty(params)
+        @statuscheck FMI.fmi3SetFloat64(wrapper.instance, wrapper.param_value_references,
+            params)
+    end
+    @statuscheck FMI.fmi3EnterInitializationMode(
+        wrapper.instance, FMI.fmi3False, zero(FMI.fmi3Float64), t, FMI.fmi3False, t)
+    if !isempty(inputs)
+        @statuscheck FMI.fmi3SetFloat64(
+            wrapper.instance, wrapper.input_value_references, inputs)
+    end
+
+    return wrapper.instance
+end
+
+"""
+    $(TYPEDSIGNATURES)
+
+Create an instance of a Model Exchange FMU. Use the existing instance in `wrapper` if
+present and create a new one otherwise. Return the instance.
+
+See `get_instance_common!` for a description of the arguments.
+"""
+function get_instance_ME!(wrapper::FMI3InstanceWrapper, inputs, params, t)
+    if wrapper.instance === nothing
+        wrapper.instance = FMI.fmi3InstantiateModelExchange!(wrapper.fmu)::FMI.FMU3Instance
+        get_instance_common!(wrapper, inputs, params, t)
+        @statuscheck FMI.fmi3ExitInitializationMode(wrapper.instance)
+        eventInfo = FMI.fmi3UpdateDiscreteStates(wrapper.instance)
+        @assert eventInfo[1] == FMI.fmi2False
+        # TODO: Support FMU events
+        @statuscheck FMI.fmi3EnterContinuousTimeMode(wrapper.instance)
+    end
+
+    return wrapper.instance
+end
+
+"""
+    $(TYPEDSIGNATURES)
+
+Create an instance of a CoSimulation FMU. Use the existing instance in `wrapper` if
+present and create a new one otherwise. Return the instance.
+
+See `get_instance_common!` for a description of the arguments.
+"""
+function get_instance_CS!(wrapper::FMI3InstanceWrapper, states, inputs, params, t)
+    if wrapper.instance === nothing
+        wrapper.instance = FMI.fmi3InstantiateCoSimulation!(
+            wrapper.fmu; eventModeUsed = false)::FMI.FMU3Instance
+        get_instance_common!(wrapper, inputs, params, t)
+        if !isempty(states)
+            @statuscheck FMI.fmi3SetFloat64(
+                wrapper.instance, wrapper.state_value_references, states)
+        end
+        @statuscheck FMI.fmi3ExitInitializationMode(wrapper.instance)
+    end
+    return wrapper.instance
+end
+
+"""
+    $(TYPEDSIGNATURES)
+"""
+function partiallyCompleteIntegratorStep(wrapper::FMI3InstanceWrapper)
+    enterEventMode = Ref(FMI.fmi3False)
+    terminateSimulation = Ref(FMI.fmi3False)
+    @statuscheck FMI.fmi3CompletedIntegratorStep!(
+        wrapper.instance, FMI.fmi3False, enterEventMode, terminateSimulation)
+    @assert enterEventMode[] == FMI.fmi3False
+    @assert terminateSimulation[] == FMI.fmi3False
+end
+
+"""
+    $(TYPEDSIGNATURES)
+"""
+function reset_instance!(wrapper::FMI3InstanceWrapper)
+    wrapper.instance === nothing && return
+    FMI.fmi3Terminate(wrapper.instance)
+    FMI.fmi3FreeInstance!(wrapper.instance)
+    wrapper.instance = nothing
+end
+
+@register_array_symbolic (fn::FMI2InstanceWrapper)(
+    states::Vector{<:Real}, inputs::Vector{<:Real}, params::Vector{<:Real}, t::Real) begin
+    size = (length(states) + length(fn.output_value_references),)
+    eltype = eltype(states)
+    ndims = 1
+end
+
+@register_array_symbolic (fn::FMI3InstanceWrapper)(
+    wrapper::FMI3InstanceWrapper, states::Vector{<:Real},
+    inputs::Vector{<:Real}, params::Vector{<:Real}, t::Real) begin
+    size = (length(states) + length(fn.output_value_references),)
+    eltype = eltype(states)
+    ndims = 1
+end
+
+"""
+    $(TYPEDSIGNATURES)
+
+Update the internal state of the ME FMU and return a vector of updated values
+for continuous state derivatives and output variables respectively. Needs to be a
+callable struct to enable symbolic registration with an inferred return size.
+"""
+function (wrapper::Union{FMI2InstanceWrapper, FMI3InstanceWrapper})(
+        states, inputs, params, t)
+    instance = get_instance_ME!(wrapper, inputs, params, t)
+
+    # TODO: Find a way to do this without allocating. We can't pass a view to these
+    # functions.
+    states_buffer = zeros(length(states))
+    outputs_buffer = zeros(length(wrapper.output_value_references))
+    # Defined in FMIBase.jl/src/eval.jl
+    # Doesn't seem to be documented, but somehow this is the only way to
+    # propagate inputs to the FMU consistently. I have no idea why.
+    instance(; x = states, u = inputs, u_refs = wrapper.input_value_references,
+        p = params, p_refs = wrapper.param_value_references, t = t)
+    # the spec requires completing the step before getting updated derivative/output values
+    partiallyCompleteIntegratorStep(wrapper)
+    instance(; dx = states_buffer, dx_refs = wrapper.derivative_value_references,
+        y = outputs_buffer, y_refs = wrapper.output_value_references)
+    return [states_buffer; outputs_buffer]
+end
+
+"""
+    $(TYPEDSIGNATURES)
+
+An affect function for use inside a `FunctionalAffect`. This should be triggered at the
+end of the solve, regardless of whether it succeeded or failed. Expects `p` to be a
+1-length array containing the index of the instance wrapper (`FMI2InstanceWrapper` or
+`FMI3InstanceWrapper`) in the parameter object.
+"""
+function fmiFinalize!(integrator, u, p, ctx)
+    wrapper_idx = p[1]
+    wrapper = integrator.ps[wrapper_idx]
+    reset_instance!(wrapper)
+end
+
+"""
+    $(TYPEDEF)
+
+A callable struct useful for initializing v2 CoSimulation FMUs. When called, updates the
+internal state of the FMU and gets updated values for output variables.
+
+# Fields
+
+$(TYPEDFIELDS)
+"""
+struct FMI2CSFunctor
+    """
+    The value references of state variables in the FMU.
+    """
+    state_value_references::Vector{FMI.fmi2ValueReference}
+    """
+    The value references of output variables in the FMU.
+    """
+    output_value_references::Vector{FMI.fmi2ValueReference}
+end
+
+function (fn::FMI2CSFunctor)(wrapper::FMI2InstanceWrapper, states, inputs, params, t)
+    states = states isa SubArray ? copy(states) : states
+    inputs = inputs isa SubArray ? copy(inputs) : inputs
+    params = params isa SubArray ? copy(params) : params
+    if wrapper.instance !== nothing
+        reset_instance!(wrapper)
+    end
+    instance = get_instance_CS!(wrapper, states, inputs, params, t)
+    if isempty(fn.output_value_references)
+        return eltype(states)[]
+    else
+        return FMI.fmi2GetReal(instance, fn.output_value_references)
+    end
+end
+
+@register_array_symbolic (fn::FMI2CSFunctor)(
+    wrapper::FMI2InstanceWrapper, states::Vector{<:Real},
+    inputs::Vector{<:Real}, params::Vector{<:Real}, t::Real) begin
+    size = (length(states) + length(fn.output_value_references),)
+    eltype = eltype(states)
+    ndims = 1
+end
+
+"""
+    $(TYPEDSIGNATURES)
+
+An affect function designed for use with `ImperativeAffect`. Should be triggered during
+callback initialization. `m` should contain the key `:states` with the value being the
+state vector if the FMU has continuous states. `m` should contain the key `:outputs` with
+the value being the output vector if the FMU has output variables. `o` should contain the
+`:inputs`, `:params`, `:t` and `:wrapper` where the latter contains the `FMI2InstanceWrapper`.
+
+Initializes the FMU. Only for use with CoSimulation FMUs.
+"""
+function fmiCSInitialize!(m, o, ctx::FMI2CSFunctor, integrator)
+    states = isdefined(m, :states) ? m.states : ()
+    inputs = o.inputs
+    params = o.params
+    t = o.t
+    wrapper = o.wrapper
+    if wrapper.instance !== nothing
+        reset_instance!(wrapper)
+    end
+
+    instance = get_instance_CS!(wrapper, states, inputs, params, t)
+    if isdefined(m, :states)
+        @statuscheck FMI.fmi2GetReal!(instance, ctx.state_value_references, m.states)
+    end
+    if isdefined(m, :outputs)
+        @statuscheck FMI.fmi2GetReal!(instance, ctx.output_value_references, m.outputs)
+    end
+
+    return m
+end
+
+"""
+    $(TYPEDSIGNATURES)
+
+An affect function designed for use with `ImperativeAffect`. Should be triggered
+periodically to communicte with the CoSimulation FMU. Has the same requirements as
+`fmiCSInitialize!` for `m` and `o`, with the addition that `o` should have a key
+`:dt` with the value being the communication step size.
+"""
+function fmiCSStep!(m, o, ctx::FMI2CSFunctor, integrator)
+    wrapper = o.wrapper
+    states = isdefined(m, :states) ? m.states : ()
+    inputs = o.inputs
+    params = o.params
+    t = o.t
+    dt = o.dt
+
+    instance = get_instance_CS!(wrapper, states, inputs, params, integrator.t)
+    if !isempty(inputs)
+        FMI.fmi2SetReal(
+            instance, wrapper.input_value_references, Csize_t(length(inputs)), inputs)
+    end
+    @statuscheck FMI.fmi2DoStep(instance, integrator.t - dt, dt, FMI.fmi2True)
+
+    if isdefined(m, :states)
+        @statuscheck FMI.fmi2GetReal!(instance, ctx.state_value_references, m.states)
+    end
+    if isdefined(m, :outputs)
+        @statuscheck FMI.fmi2GetReal!(instance, ctx.output_value_references, m.outputs)
+    end
+
+    return m
+end
+
+"""
+    $(TYPEDEF)
+
+A callable struct useful for initializing v3 CoSimulation FMUs. When called, updates the
+internal state of the FMU and gets updated values for output variables.
+
+# Fields
+
+$(TYPEDFIELDS)
+"""
+struct FMI3CSFunctor
+    """
+    The value references of state variables in the FMU.
+    """
+    state_value_references::Vector{FMI.fmi3ValueReference}
+    """
+    The value references of output variables in the FMU.
+    """
+    output_value_references::Vector{FMI.fmi3ValueReference}
+end
+
+function (fn::FMI3CSFunctor)(wrapper::FMI3InstanceWrapper, states, inputs, params, t)
+    states = states isa SubArray ? copy(states) : states
+    inputs = inputs isa SubArray ? copy(inputs) : inputs
+    params = params isa SubArray ? copy(params) : params
+    instance = get_instance_CS!(wrapper, states, inputs, params, t)
+
+    if isempty(fn.output_value_references)
+        return eltype(states)[]
+    else
+        return FMI.fmi3GetFloat64(instance, fn.output_value_references)
+    end
+end
+
+@register_array_symbolic (fn::FMI3CSFunctor)(
+    wrapper::FMI3InstanceWrapper, states::Vector{<:Real},
+    inputs::Vector{<:Real}, params::Vector{<:Real}, t::Real) begin
+    size = (length(states) + length(fn.output_value_references),)
+    eltype = eltype(states)
+    ndims = 1
+end
+
+"""
+    $(TYPEDSIGNATURES)
+"""
+function fmiCSInitialize!(m, o, ctx::FMI3CSFunctor, integrator)
+    states = isdefined(m, :states) ? m.states : ()
+    inputs = o.inputs
+    params = o.params
+    t = o.t
+    wrapper = o.wrapper
+    if wrapper.instance !== nothing
+        reset_instance!(wrapper)
+    end
+    instance = get_instance_CS!(wrapper, states, inputs, params, t)
+    if isdefined(m, :states)
+        @statuscheck FMI.fmi3GetFloat64!(instance, ctx.state_value_references, m.states)
+    end
+    if isdefined(m, :outputs)
+        @statuscheck FMI.fmi3GetFloat64!(instance, ctx.output_value_references, m.outputs)
+    end
+
+    return m
+end
+
+"""
+    $(TYPEDSIGNATURES)
+"""
+function fmiCSStep!(m, o, ctx::FMI3CSFunctor, integrator)
+    wrapper = o.wrapper
+    states = isdefined(m, :states) ? m.states : ()
+    inputs = o.inputs
+    params = o.params
+    t = o.t
+    dt = o.dt
+
+    instance = get_instance_CS!(wrapper, states, inputs, params, integrator.t)
+    if !isempty(inputs)
+        FMI.fmi3SetFloat64(instance, wrapper.input_value_references, inputs)
+    end
+    eventEncountered = Ref(FMI.fmi3False)
+    terminateSimulation = Ref(FMI.fmi3False)
+    earlyReturn = Ref(FMI.fmi3False)
+    lastSuccessfulTime = Ref(zero(FMI.fmi3Float64))
+    @statuscheck FMI.fmi3DoStep!(
+        instance, integrator.t - dt, dt, FMI.fmi3True, eventEncountered,
+        terminateSimulation, earlyReturn, lastSuccessfulTime)
+    @assert eventEncountered[] == FMI.fmi3False
+    @assert terminateSimulation[] == FMI.fmi3False
+    @assert earlyReturn[] == FMI.fmi3False
+
+    if isdefined(m, :states)
+        @statuscheck FMI.fmi3GetFloat64!(instance, ctx.state_value_references, m.states)
+    end
+    if isdefined(m, :outputs)
+        @statuscheck FMI.fmi3GetFloat64!(instance, ctx.output_value_references, m.outputs)
+    end
+
+    return m
+end
+
+end # module
diff --git a/src/ModelingToolkit.jl b/src/ModelingToolkit.jl
index 743dc8ce76..2710d7d1e4 100644
--- a/src/ModelingToolkit.jl
+++ b/src/ModelingToolkit.jl
@@ -297,4 +297,6 @@ export HomotopyContinuationProblem
 export AnalysisPoint, get_sensitivity_function, get_comp_sensitivity_function,
        get_looptransfer_function, get_sensitivity, get_comp_sensitivity, get_looptransfer,
        open_loop
+function FMIComponent end
+
 end # module
diff --git a/src/clock.jl b/src/clock.jl
index 86b7296a6d..2cd00a24bc 100644
--- a/src/clock.jl
+++ b/src/clock.jl
@@ -29,7 +29,7 @@ true if `x` contains only continuous-domain signals.
 See also [`has_continuous_domain`](@ref)
 """
 function is_continuous_domain(x)
-    issym(x) && return getmetadata(x, VariableTimeDomain, false) == Continuous
+    issym(x) && return getmetadata(x, VariableTimeDomain, false) == Continuous()
     !has_discrete_domain(x) && has_continuous_domain(x)
 end
 
@@ -58,8 +58,8 @@ has_time_domain(x::Num) = has_time_domain(value(x))
 has_time_domain(x) = false
 
 for op in [Differential]
-    @eval input_timedomain(::$op, arg = nothing) = Continuous
-    @eval output_timedomain(::$op, arg = nothing) = Continuous
+    @eval input_timedomain(::$op, arg = nothing) = Continuous()
+    @eval output_timedomain(::$op, arg = nothing) = Continuous()
 end
 
 """
diff --git a/src/discretedomain.jl b/src/discretedomain.jl
index 95abe02a7a..c7f90007c5 100644
--- a/src/discretedomain.jl
+++ b/src/discretedomain.jl
@@ -247,7 +247,7 @@ function output_timedomain(s::Shift, arg = nothing)
     InferredDiscrete
 end
 
-input_timedomain(::Sample, _ = nothing) = Continuous
+input_timedomain(::Sample, _ = nothing) = Continuous()
 output_timedomain(s::Sample, _ = nothing) = s.clock
 
 function input_timedomain(h::Hold, arg = nothing)
@@ -256,7 +256,7 @@ function input_timedomain(h::Hold, arg = nothing)
     end
     InferredDiscrete # the Hold accepts any discrete
 end
-output_timedomain(::Hold, _ = nothing) = Continuous
+output_timedomain(::Hold, _ = nothing) = Continuous()
 
 sampletime(op::Sample, _ = nothing) = sampletime(op.clock)
 sampletime(op::ShiftIndex, _ = nothing) = sampletime(op.clock)
diff --git a/src/parameters.jl b/src/parameters.jl
index ced7767718..91121b7cbb 100644
--- a/src/parameters.jl
+++ b/src/parameters.jl
@@ -33,7 +33,12 @@ end
 
 function getcalledparameter(x)
     x = unwrap(x)
-    return getmetadata(x, CallWithParent)
+    # `parent` is a `CallWithMetadata` with the correct metadata,
+    # but no namespacing. `operation(x)` has the correct namespacing,
+    # but is not a `CallWithMetadata` and doesn't have any metadata.
+    # This approach combines both.
+    parent = getmetadata(x, CallWithParent)
+    return CallWithMetadata(operation(x), metadata(parent))
 end
 
 """
diff --git a/src/structural_transformation/symbolics_tearing.jl b/src/structural_transformation/symbolics_tearing.jl
index 8161a9572b..d95951c41e 100644
--- a/src/structural_transformation/symbolics_tearing.jl
+++ b/src/structural_transformation/symbolics_tearing.jl
@@ -593,7 +593,7 @@ function tearing_reassemble(state::TearingState, var_eq_matching,
     @set! sys.unknowns = unknowns
 
     obs, subeqs, deps = cse_and_array_hacks(
-        obs, subeqs, unknowns, neweqs; cse = cse_hack, array = array_hack)
+        sys, obs, subeqs, unknowns, neweqs; cse = cse_hack, array = array_hack)
 
     @set! sys.eqs = neweqs
     @set! sys.observed = obs
@@ -627,7 +627,7 @@ if all `p[i]` are present and the unscalarized form is used in any equation (obs
 not) we first count the number of times the scalarized form of each observed variable
 occurs in observed equations (and unknowns if it's split).
 """
-function cse_and_array_hacks(obs, subeqs, unknowns, neweqs; cse = true, array = true)
+function cse_and_array_hacks(sys, obs, subeqs, unknowns, neweqs; cse = true, array = true)
     # HACK 1
     # mapping of rhs to temporary CSE variable
     # `f(...) => tmpvar` in above example
@@ -696,6 +696,7 @@ function cse_and_array_hacks(obs, subeqs, unknowns, neweqs; cse = true, array =
                     tempvar; T = Symbolics.symtype(rhs_arr)))
                 tempvar = setmetadata(
                     tempvar, Symbolics.ArrayShapeCtx, Symbolics.shape(rhs_arr))
+                vars!(all_vars, rhs_arr)
                 tempeq = tempvar ~ rhs_arr
                 rhs_to_tempvar[rhs_arr] = tempvar
                 push!(obs, tempeq)
@@ -718,12 +719,16 @@ function cse_and_array_hacks(obs, subeqs, unknowns, neweqs; cse = true, array =
         Symbolics.shape(sym) != Symbolics.Unknown() || continue
         arg1 = arguments(sym)[1]
         cnt = get(arr_obs_occurrences, arg1, 0)
-        cnt == 0 && continue
         arr_obs_occurrences[arg1] = cnt + 1
     end
     for eq in neweqs
         vars!(all_vars, eq.rhs)
     end
+
+    # also count unscalarized variables used in callbacks
+    for ev in Iterators.flatten((continuous_events(sys), discrete_events(sys)))
+        vars!(all_vars, ev)
+    end
     obs_arr_eqs = Equation[]
     for (arrvar, cnt) in arr_obs_occurrences
         cnt == length(arrvar) || continue
@@ -737,7 +742,9 @@ function cse_and_array_hacks(obs, subeqs, unknowns, neweqs; cse = true, array =
         # try to `create_array(OffsetArray{...}, ...)` which errors.
         # `term(Origin(firstind), scal)` doesn't retain the `symtype` and `size`
         # of `scal`.
-        push!(obs_arr_eqs, arrvar ~ change_origin(Origin(firstind), scal))
+        rhs = scal
+        rhs = change_origin(firstind, rhs)
+        push!(obs_arr_eqs, arrvar ~ rhs)
     end
     append!(obs, obs_arr_eqs)
     append!(subeqs, obs_arr_eqs)
@@ -764,10 +771,13 @@ getindex_wrapper(x, i) = x[i...]
 
 # PART OF HACK 2
 function change_origin(origin, arr)
-    return origin(arr)
+    if all(isone, Tuple(origin))
+        return arr
+    end
+    return Origin(origin)(arr)
 end
 
-@register_array_symbolic change_origin(origin::Origin, arr::AbstractArray) begin
+@register_array_symbolic change_origin(origin::Any, arr::AbstractArray) begin
     size = size(arr)
     eltype = eltype(arr)
     ndims = ndims(arr)
diff --git a/src/systems/callbacks.jl b/src/systems/callbacks.jl
index eaf31a9c5f..492b5dac8c 100644
--- a/src/systems/callbacks.jl
+++ b/src/systems/callbacks.jl
@@ -77,6 +77,13 @@ function has_functional_affect(cb)
     (affects(cb) isa FunctionalAffect || affects(cb) isa ImperativeAffect)
 end
 
+function vars!(vars, aff::FunctionalAffect; op = Differential)
+    for var in Iterators.flatten((unknowns(aff), parameters(aff), discretes(aff)))
+        vars!(vars, var)
+    end
+    return vars
+end
+
 #################################### continuous events #####################################
 
 const NULL_AFFECT = Equation[]
@@ -333,6 +340,22 @@ function continuous_events(sys::AbstractSystem)
     filter(!isempty, cbs)
 end
 
+function vars!(vars, cb::SymbolicContinuousCallback; op = Differential)
+    for eq in equations(cb)
+        vars!(vars, eq; op)
+    end
+    for aff in (affects(cb), affect_negs(cb), initialize_affects(cb), finalize_affects(cb))
+        if aff isa Vector{Equation}
+            for eq in aff
+                vars!(vars, eq; op)
+            end
+        elseif aff !== nothing
+            vars!(vars, aff; op)
+        end
+    end
+    return vars
+end
+
 #################################### discrete events #####################################
 
 struct SymbolicDiscreteCallback
@@ -469,6 +492,28 @@ function discrete_events(sys::AbstractSystem)
     cbs
 end
 
+function vars!(vars, cb::SymbolicDiscreteCallback; op = Differential)
+    if symbolic_type(cb.condition) == NotSymbolic
+        if cb.condition isa AbstractArray
+            for eq in cb.condition
+                vars!(vars, eq; op)
+            end
+        end
+    else
+        vars!(vars, cb.condition; op)
+    end
+    for aff in (cb.affects, cb.initialize, cb.finalize)
+        if aff isa Vector{Equation}
+            for eq in aff
+                vars!(vars, eq; op)
+            end
+        elseif aff !== nothing
+            vars!(vars, aff; op)
+        end
+    end
+    return vars
+end
+
 ################################# compilation functions ####################################
 
 # handles ensuring that affect! functions work with integrator arguments
diff --git a/src/systems/clock_inference.jl b/src/systems/clock_inference.jl
index a92b2aa67c..611f8e2fae 100644
--- a/src/systems/clock_inference.jl
+++ b/src/systems/clock_inference.jl
@@ -8,8 +8,8 @@ end
 function ClockInference(ts::TransformationState)
     @unpack structure = ts
     @unpack graph = structure
-    eq_domain = TimeDomain[Continuous for _ in 1:nsrcs(graph)]
-    var_domain = TimeDomain[Continuous for _ in 1:ndsts(graph)]
+    eq_domain = TimeDomain[Continuous() for _ in 1:nsrcs(graph)]
+    var_domain = TimeDomain[Continuous() for _ in 1:ndsts(graph)]
     inferred = BitSet()
     for (i, v) in enumerate(get_fullvars(ts))
         d = get_time_domain(ts, v)
@@ -151,7 +151,7 @@ function split_system(ci::ClockInference{S}) where {S}
             get!(clock_to_id, d) do
                 cid = (cid_counter[] += 1)
                 push!(id_to_clock, d)
-                if d == Continuous
+                if d == Continuous()
                     continuous_id[] = cid
                 end
                 cid
diff --git a/src/systems/diffeqs/abstractodesystem.jl b/src/systems/diffeqs/abstractodesystem.jl
index ac28d41dd4..6571aed9cc 100644
--- a/src/systems/diffeqs/abstractodesystem.jl
+++ b/src/systems/diffeqs/abstractodesystem.jl
@@ -1391,14 +1391,16 @@ function InitializationProblem{iip, specialize}(sys::AbstractSystem,
     fullmap = merge(u0map, parammap)
     u0T = Union{}
     for sym in unknowns(isys)
-        haskey(fullmap, sym) || continue
-        symbolic_type(fullmap[sym]) == NotSymbolic() || continue
-        u0T = promote_type(u0T, typeof(fullmap[sym]))
+        val = fixpoint_sub(sym, fullmap)
+        symbolic_type(val) == NotSymbolic() || continue
+        u0T = promote_type(u0T, typeof(val))
     end
     for eq in observed(isys)
-        haskey(fullmap, eq.lhs) || continue
-        symbolic_type(fullmap[eq.lhs]) == NotSymbolic() || continue
-        u0T = promote_type(u0T, typeof(fullmap[eq.lhs]))
+        # ignore HACK-ed observed equations
+        symbolic_type(eq.lhs) == ArraySymbolic() && continue
+        val = fixpoint_sub(eq.lhs, fullmap)
+        symbolic_type(val) == NotSymbolic() || continue
+        u0T = promote_type(u0T, typeof(val))
     end
     if u0T != Union{}
         u0T = eltype(u0T)
diff --git a/src/systems/diffeqs/odesystem.jl b/src/systems/diffeqs/odesystem.jl
index d9e126baba..41bcaf1d01 100644
--- a/src/systems/diffeqs/odesystem.jl
+++ b/src/systems/diffeqs/odesystem.jl
@@ -373,6 +373,8 @@ function Base.:(==)(sys1::ODESystem, sys2::ODESystem)
         _eq_unordered(get_eqs(sys1), get_eqs(sys2)) &&
         _eq_unordered(get_unknowns(sys1), get_unknowns(sys2)) &&
         _eq_unordered(get_ps(sys1), get_ps(sys2)) &&
+        _eq_unordered(continuous_events(sys1), continuous_events(sys2)) &&
+        _eq_unordered(discrete_events(sys1), discrete_events(sys2)) && 
         all(s1 == s2 for (s1, s2) in zip(get_systems(sys1), get_systems(sys2)))
 end
 
@@ -685,13 +687,22 @@ function populate_delays(delays::Set, obsexprs, histfn, sys, sym)
 end
 
 function _eq_unordered(a, b)
+    # a and b may be multidimensional
+    # e.g. comparing noiseeqs of SDESystem
+    a = vec(a)
+    b = vec(b)
     length(a) === length(b) || return false
     n = length(a)
     idxs = Set(1:n)
     for x in a
         idx = findfirst(isequal(x), b)
+        # loop since there might be multiple identical entries in a/b
+        # and while we might have already matched the first there could
+        # be a second that is equal to x
+        while idx !== nothing && !(idx in idxs)
+            idx = findnext(isequal(x), b, idx + 1)
+        end
         idx === nothing && return false
-        idx ∈ idxs || return false
         delete!(idxs, idx)
     end
     return true
diff --git a/src/systems/diffeqs/sdesystem.jl b/src/systems/diffeqs/sdesystem.jl
index fb2eb1d804..df950ee01c 100644
--- a/src/systems/diffeqs/sdesystem.jl
+++ b/src/systems/diffeqs/sdesystem.jl
@@ -279,11 +279,13 @@ function Base.:(==)(sys1::SDESystem, sys2::SDESystem)
     iv2 = get_iv(sys2)
     isequal(iv1, iv2) &&
         isequal(nameof(sys1), nameof(sys2)) &&
-        isequal(get_eqs(sys1), get_eqs(sys2)) &&
-        isequal(get_noiseeqs(sys1), get_noiseeqs(sys2)) &&
+        _eq_unordered(get_eqs(sys1), get_eqs(sys2)) &&
+        _eq_unordered(get_noiseeqs(sys1), get_noiseeqs(sys2)) &&
         isequal(get_is_scalar_noise(sys1), get_is_scalar_noise(sys2)) &&
         _eq_unordered(get_unknowns(sys1), get_unknowns(sys2)) &&
         _eq_unordered(get_ps(sys1), get_ps(sys2)) &&
+        _eq_unordered(continuous_events(sys1), continuous_events(sys2)) &&
+        _eq_unordered(discrete_events(sys1), discrete_events(sys2)) && 
         all(s1 == s2 for (s1, s2) in zip(get_systems(sys1), get_systems(sys2)))
 end
 
diff --git a/src/systems/imperative_affect.jl b/src/systems/imperative_affect.jl
index 2f489913ad..aee95fd051 100644
--- a/src/systems/imperative_affect.jl
+++ b/src/systems/imperative_affect.jl
@@ -38,7 +38,7 @@ in the returned tuple, in which case the associated field will not be updated.
     skip_checks::Bool
 end
 
-function ImperativeAffect(f::Function;
+function ImperativeAffect(f;
         observed::NamedTuple = NamedTuple{()}(()),
         modified::NamedTuple = NamedTuple{()}(()),
         ctx = nothing,
@@ -48,18 +48,18 @@ function ImperativeAffect(f::Function;
         collect(values(modified)), collect(keys(modified)),
         ctx, skip_checks)
 end
-function ImperativeAffect(f::Function, modified::NamedTuple;
+function ImperativeAffect(f, modified::NamedTuple;
         observed::NamedTuple = NamedTuple{()}(()), ctx = nothing, skip_checks = false)
     ImperativeAffect(
         f, observed = observed, modified = modified, ctx = ctx, skip_checks = skip_checks)
 end
 function ImperativeAffect(
-        f::Function, modified::NamedTuple, observed::NamedTuple; ctx = nothing, skip_checks = false)
+        f, modified::NamedTuple, observed::NamedTuple; ctx = nothing, skip_checks = false)
     ImperativeAffect(
         f, observed = observed, modified = modified, ctx = ctx, skip_checks = skip_checks)
 end
 function ImperativeAffect(
-        f::Function, modified::NamedTuple, observed::NamedTuple, ctx; skip_checks = false)
+        f, modified::NamedTuple, observed::NamedTuple, ctx; skip_checks = false)
     ImperativeAffect(
         f, observed = observed, modified = modified, ctx = ctx, skip_checks = skip_checks)
 end
@@ -216,3 +216,20 @@ function compile_user_affect(affect::ImperativeAffect, cb, sys, dvs, ps; kwargs.
 end
 
 scalarize_affects(affects::ImperativeAffect) = affects
+
+function vars!(vars, aff::ImperativeAffect; op = Differential)
+    for var in Iterators.flatten((observed(aff), modified(aff)))
+        if symbolic_type(var) == NotSymbolic()
+            if var isa AbstractArray
+                for v in var
+                    v = unwrap(v)
+                    vars!(vars, v)
+                end
+            end
+        else
+            var = unwrap(var)
+            vars!(vars, var)
+        end
+    end
+    return vars
+end
diff --git a/src/systems/index_cache.jl b/src/systems/index_cache.jl
index 623623da42..e4ca0b6b48 100644
--- a/src/systems/index_cache.jl
+++ b/src/systems/index_cache.jl
@@ -299,6 +299,9 @@ function IndexCache(sys::AbstractSystem)
                 for v in vs
                     if (idx = get(disc_idxs, v, nothing)) !== nothing
                         push!(timeseries, idx.clock_idx)
+                    elseif iscall(v) && operation(v) === getindex &&
+                           (idx = get(disc_idxs, arguments(v)[1], nothing)) !== nothing
+                        push!(timeseries, idx.clock_idx)
                     elseif haskey(observed_syms_to_timeseries, v)
                         union!(timeseries, observed_syms_to_timeseries[v])
                     elseif haskey(dependent_pars_to_timeseries, v)
diff --git a/src/systems/nonlinear/initializesystem.jl b/src/systems/nonlinear/initializesystem.jl
index a7f98d79b1..296078bf43 100644
--- a/src/systems/nonlinear/initializesystem.jl
+++ b/src/systems/nonlinear/initializesystem.jl
@@ -193,6 +193,15 @@ function generate_initializesystem(sys::AbstractSystem;
         append!(eqs_ics, trueobs)
     end
 
+    # even if `p => tovar(p)` is in `paramsubs`, `isparameter(p[1]) === true` after substitution
+    # so add scalarized versions as well
+    for k in collect(keys(paramsubs))
+        symbolic_type(k) == ArraySymbolic() || continue
+        for i in eachindex(k)
+            paramsubs[k[i]] = paramsubs[k][i]
+        end
+    end
+
     eqs_ics = Symbolics.substitute.(eqs_ics, (paramsubs,))
     if is_time_dependent(sys)
         vars = [vars; collect(values(paramsubs))]
@@ -223,15 +232,31 @@ struct ReconstructInitializeprob
 end
 
 function ReconstructInitializeprob(srcsys::AbstractSystem, dstsys::AbstractSystem)
-    syms = [unknowns(dstsys);
-            reduce(vcat, reorder_parameters(dstsys, parameters(dstsys)); init = [])]
+    syms = reduce(vcat, reorder_parameters(dstsys, parameters(dstsys)); init = [])
     getter = getu(srcsys, syms)
-    setter = setsym_oop(dstsys, syms)
+    setter = setp_oop(dstsys, syms)
     return ReconstructInitializeprob(getter, setter)
 end
 
 function (rip::ReconstructInitializeprob)(srcvalp, dstvalp)
-    rip.setter(dstvalp, rip.getter(srcvalp))
+    newp = rip.setter(dstvalp, rip.getter(srcvalp))
+    if state_values(dstvalp) === nothing
+        return nothing, newp
+    end
+    T = eltype(state_values(srcvalp))
+    if parameter_values(dstvalp) isa MTKParameters
+        if !isempty(newp.tunable)
+            T = promote_type(eltype(newp.tunable), T)
+        end
+    elseif !isempty(newp)
+        T = promote_type(eltype(newp), T)
+    end
+    if T == eltype(state_values(dstvalp))
+        u0 = state_values(dstvalp)
+    else
+        u0 = T.(state_values(dstvalp))
+    end
+    return u0, newp
 end
 
 struct InitializationSystemMetadata
diff --git a/src/systems/nonlinear/nonlinearsystem.jl b/src/systems/nonlinear/nonlinearsystem.jl
index 6b0b0cc759..89663f7dbc 100644
--- a/src/systems/nonlinear/nonlinearsystem.jl
+++ b/src/systems/nonlinear/nonlinearsystem.jl
@@ -710,7 +710,7 @@ function SciMLBase.SCCNonlinearProblem{iip}(sys::NonlinearSystem, u0map,
         # precomputed subexpressions should not contain `banned_vars`
         banned_vars = Set{Any}(vcat(_dvs, getproperty.(_obs, (:lhs,))))
         filter!(banned_vars) do var
-            symbolic_type(var) != ArraySymbolic() || all(x -> var[i] in banned_vars, eachindex(var))
+            symbolic_type(var) != ArraySymbolic() || all(j -> var[j] in banned_vars, eachindex(var))
         end
         state = Dict()
         for i in eachindex(_obs)
@@ -772,6 +772,7 @@ function SciMLBase.SCCNonlinearProblem{iip}(sys::NonlinearSystem, u0map,
         _obs = scc_obs[i]
         cachevars = scc_cachevars[i]
         cacheexprs = scc_cacheexprs[i]
+        _prevobsidxs = vcat(_prevobsidxs, observed_equations_used_by(sys, reduce(vcat, values(cacheexprs); init = [])))
 
         if isempty(cachevars)
             push!(explicitfuns, Returns(nothing))
diff --git a/src/systems/problem_utils.jl b/src/systems/problem_utils.jl
index a8d443349b..81be9c3d9e 100644
--- a/src/systems/problem_utils.jl
+++ b/src/systems/problem_utils.jl
@@ -88,6 +88,18 @@ function symbols_to_symbolics!(sys::AbstractSystem, varmap::AbstractDict)
     end
 end
 
+"""
+    $(TYPEDSIGNATURES)
+
+Utility function to get the value `val` corresponding to key `var` in `varmap`, and
+return `getindex(val, idx)` if it exists or `nothing` otherwise.
+"""
+function get_and_getindex(varmap, var, idx)
+    val = get(varmap, var, nothing)
+    val === nothing && return nothing
+    return val[idx]
+end
+
 """
     $(TYPEDSIGNATURES)
 
@@ -115,8 +127,9 @@ function add_fallbacks!(
             val = map(eachindex(var)) do idx
                 # @something is lazy and saves from writing a massive if-elseif-else
                 @something(get(varmap, var[idx], nothing),
-                    get(varmap, ttvar[idx], nothing), get(fallbacks, var, nothing)[idx],
-                    get(fallbacks, ttvar, nothing)[idx], get(fallbacks, var[idx], nothing),
+                    get(varmap, ttvar[idx], nothing), get_and_getindex(fallbacks, var, idx),
+                    get_and_getindex(fallbacks, ttvar, idx), get(
+                        fallbacks, var[idx], nothing),
                     get(fallbacks, ttvar[idx], nothing), Some(nothing))
             end
             # only push the missing entries
@@ -578,6 +591,10 @@ function maybe_build_initialization_problem(
             p = unwrap(p)
             stype = symtype(p)
             op[p] = get_temporary_value(p)
+            if iscall(p) && operation(p) === getindex
+                arrp = arguments(p)[1]
+                op[arrp] = collect(arrp)
+            end
         end
 
         if is_time_dependent(sys)
diff --git a/src/systems/systems.jl b/src/systems/systems.jl
index 04c50bc766..9c8c272c5c 100644
--- a/src/systems/systems.jl
+++ b/src/systems/systems.jl
@@ -57,6 +57,7 @@ function structural_simplify(
         ks = collect(keys(defs))  # take copy to avoid mutating defs while iterating.
         for k in ks
             if Symbolics.isarraysymbolic(k) && Symbolics.shape(k) !== Symbolics.Unknown()
+                defs[k] === missing && continue
                 for i in eachindex(k)
                     defs[k[i]] = defs[k][i]
                 end
diff --git a/src/systems/systemstructure.jl b/src/systems/systemstructure.jl
index 1c0ad8b8ae..1bdc11f06a 100644
--- a/src/systems/systemstructure.jl
+++ b/src/systems/systemstructure.jl
@@ -662,7 +662,7 @@ function structural_simplify!(state::TearingState, io = nothing; simplify = fals
                 Dict(v => 0.0 for v in Iterators.flatten(inputs)))
         end
         ps = [sym isa CallWithMetadata ? sym :
-              setmetadata(sym, VariableTimeDomain, get(time_domains, sym, Continuous))
+              setmetadata(sym, VariableTimeDomain, get(time_domains, sym, Continuous()))
               for sym in get_ps(sys)]
         @set! sys.ps = ps
     else
diff --git a/src/utils.jl b/src/utils.jl
index 5e4f0b52d2..b6544972b9 100644
--- a/src/utils.jl
+++ b/src/utils.jl
@@ -403,6 +403,9 @@ function vars!(vars, eq::Equation; op = Differential)
 end
 function vars!(vars, O; op = Differential)
     if isvariable(O)
+        if iscall(O) && operation(O) === getindex && iscalledparameter(first(arguments(O)))
+            O = first(arguments(O))
+        end
         if iscalledparameter(O)
             f = getcalledparameter(O)
             push!(vars, f)
diff --git a/test/clock.jl b/test/clock.jl
index 91aaa8248e..67f469e1a6 100644
--- a/test/clock.jl
+++ b/test/clock.jl
@@ -77,19 +77,19 @@ k = ShiftIndex(d)
 
 d = Clock(dt)
 # Note that TearingState reorders the equations
-@test eqmap[1] == Continuous
+@test eqmap[1] == Continuous()
 @test eqmap[2] == d
 @test eqmap[3] == d
 @test eqmap[4] == d
-@test eqmap[5] == Continuous
-@test eqmap[6] == Continuous
+@test eqmap[5] == Continuous()
+@test eqmap[6] == Continuous()
 
 @test varmap[yd] == d
 @test varmap[ud] == d
 @test varmap[r] == d
-@test varmap[x] == Continuous
-@test varmap[y] == Continuous
-@test varmap[u] == Continuous
+@test varmap[x] == Continuous()
+@test varmap[y] == Continuous()
+@test varmap[u] == Continuous()
 
 @info "Testing shift normalization"
 dt = 0.1
@@ -192,10 +192,10 @@ eqs = [yd ~ Sample(dt)(y)
     @test varmap[ud1] == d
     @test varmap[yd2] == d2
     @test varmap[ud2] == d2
-    @test varmap[r] == Continuous
-    @test varmap[x] == Continuous
-    @test varmap[y] == Continuous
-    @test varmap[u] == Continuous
+    @test varmap[r] == Continuous()
+    @test varmap[x] == Continuous()
+    @test varmap[y] == Continuous()
+    @test varmap[u] == Continuous()
 
     @info "test composed systems"
 
@@ -241,14 +241,14 @@ eqs = [yd ~ Sample(dt)(y)
     ci, varmap = infer_clocks(cl)
 
     @test varmap[f.x] == Clock(0.5)
-    @test varmap[p.x] == Continuous
-    @test varmap[p.y] == Continuous
+    @test varmap[p.x] == Continuous()
+    @test varmap[p.y] == Continuous()
     @test varmap[c.ud] == Clock(0.5)
     @test varmap[c.yd] == Clock(0.5)
-    @test varmap[c.y] == Continuous
+    @test varmap[c.y] == Continuous()
     @test varmap[f.y] == Clock(0.5)
     @test varmap[f.u] == Clock(0.5)
-    @test varmap[p.u] == Continuous
+    @test varmap[p.u] == Continuous()
     @test varmap[c.r] == Clock(0.5)
 
     ## Multiple clock rates
@@ -474,7 +474,7 @@ eqs = [yd ~ Sample(dt)(y)
 
     ## Test continuous clock
 
-    c = ModelingToolkit.SolverStepClock
+    c = ModelingToolkit.SolverStepClock()
     k = ShiftIndex(c)
 
     @mtkmodel CounterSys begin
diff --git a/test/fmi/Project.toml b/test/fmi/Project.toml
new file mode 100644
index 0000000000..59bd3c90da
--- /dev/null
+++ b/test/fmi/Project.toml
@@ -0,0 +1,8 @@
+[deps]
+FMI = "14a09403-18e3-468f-ad8a-74f8dda2d9ac"
+FMIZoo = "724179cf-c260-40a9-bd27-cccc6fe2f195"
+ModelingToolkit = "961ee093-0014-501f-94e3-6117800e7a78"
+OrdinaryDiffEq = "1dea7af3-3e70-54e6-95c3-0bf5283fa5ed"
+
+[compat]
+FMI = "0.14"
diff --git a/test/fmi/fmi.jl b/test/fmi/fmi.jl
new file mode 100644
index 0000000000..bf5b60459f
--- /dev/null
+++ b/test/fmi/fmi.jl
@@ -0,0 +1,309 @@
+using ModelingToolkit, FMI, FMIZoo, OrdinaryDiffEq, NonlinearSolve, SciMLBase
+using ModelingToolkit: t_nounits as t, D_nounits as D
+import ModelingToolkit as MTK
+
+const FMU_DIR = joinpath(@__DIR__, "fmus")
+
+@testset "Standalone pendulum model" begin
+    fmu = loadFMU("SpringPendulum1D", "Dymola", "2022x"; type = :ME)
+    truesol = FMI.simulate(
+        fmu, (0.0, 8.0); saveat = 0.0:0.1:8.0, recordValues = ["mass.s", "mass.v"])
+
+    function test_no_inputs_outputs(sys)
+        for var in unknowns(sys)
+            @test !MTK.isinput(var)
+            @test !MTK.isoutput(var)
+        end
+    end
+    @testset "v2, ME" begin
+        fmu = loadFMU("SpringPendulum1D", "Dymola", "2022x"; type = :ME)
+        @mtkbuild sys = MTK.FMIComponent(Val(2); fmu, type = :ME)
+        test_no_inputs_outputs(sys)
+        prob = ODEProblem{true, SciMLBase.FullSpecialize}(
+            sys, [sys.mass__s => 0.5, sys.mass__v => 0.0], (0.0, 8.0))
+        sol = solve(prob, Tsit5(); reltol = 1e-8, abstol = 1e-8)
+        @test SciMLBase.successful_retcode(sol)
+
+        @test sol(0.0:0.1:8.0;
+            idxs = [sys.mass__s, sys.mass__v]).u≈collect.(truesol.values.saveval) atol=1e-4
+        # repeated solve works
+        @test_nowarn solve(prob, Tsit5())
+    end
+    @testset "v2, CS" begin
+        fmu = loadFMU("SpringPendulum1D", "Dymola", "2022x"; type = :CS)
+        @named inner = MTK.FMIComponent(
+            Val(2); fmu, communication_step_size = 0.001, type = :CS)
+        @variables x(t) = 1.0
+        @mtkbuild sys = ODESystem([D(x) ~ x], t; systems = [inner])
+        test_no_inputs_outputs(sys)
+
+        prob = ODEProblem{true, SciMLBase.FullSpecialize}(
+            sys, [sys.inner.mass__s => 0.5, sys.inner.mass__v => 0.0], (0.0, 8.0))
+        sol = solve(prob, Tsit5(); reltol = 1e-8, abstol = 1e-8)
+        @test SciMLBase.successful_retcode(sol)
+
+        @test sol(0.0:0.1:8.0;
+            idxs = [sys.inner.mass__s, sys.inner.mass__v]).u≈collect.(truesol.values.saveval) rtol=1e-2
+    end
+
+    fmu = loadFMU("SpringPendulum1D", "Dymola", "2023x", "3.0"; type = :ME)
+    truesol = FMI.simulate(
+        fmu, (0.0, 8.0); saveat = 0.0:0.1:8.0, recordValues = ["mass.s", "mass.v"])
+    @testset "v3, ME" begin
+        fmu = loadFMU("SpringPendulum1D", "Dymola", "2023x", "3.0"; type = :ME)
+        @mtkbuild sys = MTK.FMIComponent(Val(3); fmu, type = :ME)
+        test_no_inputs_outputs(sys)
+        prob = ODEProblem{true, SciMLBase.FullSpecialize}(
+            sys, [sys.mass__s => 0.5, sys.mass__v => 0.0], (0.0, 8.0))
+        sol = solve(prob, Tsit5(); reltol = 1e-8, abstol = 1e-8)
+        @test SciMLBase.successful_retcode(sol)
+
+        @test sol(0.0:0.1:8.0;
+            idxs = [sys.mass__s, sys.mass__v]).u≈collect.(truesol.values.saveval) atol=1e-4
+        # repeated solve works
+        @test_nowarn solve(prob, Tsit5())
+    end
+    @testset "v3, CS" begin
+        fmu = loadFMU("SpringPendulum1D", "Dymola", "2023x", "3.0"; type = :CS)
+        @named inner = MTK.FMIComponent(
+            Val(3); fmu, communication_step_size = 0.001, type = :CS)
+        @variables x(t) = 1.0
+        @mtkbuild sys = ODESystem([D(x) ~ x], t; systems = [inner])
+        test_no_inputs_outputs(sys)
+
+        prob = ODEProblem{true, SciMLBase.FullSpecialize}(
+            sys, [sys.inner.mass__s => 0.5, sys.inner.mass__v => 0.0], (0.0, 8.0))
+        sol = solve(prob, Tsit5(); reltol = 1e-8, abstol = 1e-8)
+        @test SciMLBase.successful_retcode(sol)
+
+        @test sol(0.0:0.1:8.0;
+            idxs = [sys.inner.mass__s, sys.inner.mass__v]).u≈collect.(truesol.values.saveval) rtol=1e-2
+    end
+end
+
+@mtkmodel SimpleAdder begin
+    @variables begin
+        a(t)
+        b(t)
+        c(t)
+        out(t)
+        out2(t)
+    end
+    @parameters begin
+        value = 1.0
+    end
+    @equations begin
+        out ~ a + b + value
+        D(c) ~ out
+        out2 ~ 2c
+    end
+end
+
+@mtkmodel StateSpace begin
+    @variables begin
+        x(t)
+        y(t)
+        u(t)
+    end
+    @parameters begin
+        A = 1.0
+        B = 1.0
+        C = 1.0
+        _D = 1.0
+    end
+    @equations begin
+        D(x) ~ A * x + B * u
+        y ~ C * x + _D * u
+    end
+end
+
+@testset "IO Model" begin
+    function build_simple_adder(adder)
+        @variables a(t) b(t) c(t) [guess = 1.0]
+        @mtkbuild sys = ODESystem(
+            [adder.a ~ a, adder.b ~ b, D(a) ~ t,
+                D(b) ~ adder.out + adder.c, c^2 ~ adder.out + adder.value],
+            t;
+            systems = [adder])
+        # c will be solved for by initialization
+        # this tests that initialization also works with FMUs
+        prob = ODEProblem(sys, [sys.adder.c => 2.0, sys.a => 1.0, sys.b => 1.0],
+            (0.0, 1.0), [sys.adder.value => 2.0])
+        return sys, prob
+    end
+
+    @named adder = SimpleAdder()
+    truesys, trueprob = build_simple_adder(adder)
+    truesol = solve(trueprob, abstol = 1e-8, reltol = 1e-8)
+    @test SciMLBase.successful_retcode(truesol)
+
+    @testset "v2, ME" begin
+        fmu = loadFMU(joinpath(FMU_DIR, "SimpleAdder.fmu"); type = :ME)
+        @named adder = MTK.FMIComponent(Val(2); fmu, type = :ME)
+        @test MTK.isinput(adder.a)
+        @test MTK.isinput(adder.b)
+        @test MTK.isoutput(adder.out)
+        @test MTK.isoutput(adder.out2)
+        @test !MTK.isinput(adder.c) && !MTK.isoutput(adder.c)
+
+        sys, prob = build_simple_adder(adder)
+        sol = solve(prob, Rodas5P(autodiff = false), abstol = 1e-8, reltol = 1e-8)
+        @test SciMLBase.successful_retcode(sol)
+
+        @test truesol(sol.t;
+            idxs = [truesys.a, truesys.b, truesys.c, truesys.adder.c]).u≈sol[[
+            sys.a, sys.b, sys.c, sys.adder.c]] rtol=1e-7
+    end
+    @testset "v2, CS" begin
+        fmu = loadFMU(joinpath(FMU_DIR, "SimpleAdder.fmu"); type = :CS)
+        @named adder = MTK.FMIComponent(
+            Val(2); fmu, type = :CS, communication_step_size = 1e-3,
+            reinitializealg = BrownFullBasicInit())
+        @test MTK.isinput(adder.a)
+        @test MTK.isinput(adder.b)
+        @test MTK.isoutput(adder.out)
+        @test MTK.isoutput(adder.out2)
+        @test !MTK.isinput(adder.c) && !MTK.isoutput(adder.c)
+
+        sys, prob = build_simple_adder(adder)
+        sol = solve(prob, Rodas5P(autodiff = false), abstol = 1e-8, reltol = 1e-8)
+        @test SciMLBase.successful_retcode(sol)
+
+        @test truesol(sol.t; idxs = [truesys.a, truesys.b, truesys.c]).u≈sol.u rtol=1e-2
+        # sys.adder.c is a discrete variable
+        @test truesol(sol.t; idxs = truesys.adder.c).u≈sol(sol.t; idxs = sys.adder.c).u rtol=1e-3
+    end
+
+    function build_sspace_model(sspace)
+        @variables u(t)=1.0 x(t)=1.0 y(t) [guess = 1.0]
+        @mtkbuild sys = ODESystem(
+            [sspace.u ~ u, D(u) ~ t, D(x) ~ sspace.x + sspace.y, y^2 ~ sspace.y + sspace.x], t;
+            systems = [sspace]
+        )
+
+        prob = ODEProblem(
+            sys, [sys.sspace.x => 1.0], (0.0, 1.0), [sys.sspace.A => 2.0]; use_scc = false)
+        return sys, prob
+    end
+
+    @named sspace = StateSpace()
+    truesys, trueprob = build_sspace_model(sspace)
+    truesol = solve(trueprob, abstol = 1e-8, reltol = 1e-8)
+    @test SciMLBase.successful_retcode(truesol)
+
+    @testset "v3, ME" begin
+        fmu = loadFMU(joinpath(FMU_DIR, "StateSpace.fmu"); type = :ME)
+        @named sspace = MTK.FMIComponent(Val(3); fmu, type = :ME)
+        @test MTK.isinput(sspace.u)
+        @test MTK.isoutput(sspace.y)
+        @test !MTK.isinput(sspace.x) && !MTK.isoutput(sspace.x)
+
+        sys, prob = build_sspace_model(sspace)
+        sol = solve(prob, Rodas5P(autodiff = false); abstol = 1e-8, reltol = 1e-8)
+        @test SciMLBase.successful_retcode(sol)
+
+        @test truesol(sol.t;
+            idxs = [truesys.u, truesys.x, truesys.y, truesys.sspace.x]).u≈sol[[
+            sys.u, sys.x, sys.y, sys.sspace.x]] rtol=1e-7
+    end
+
+    @testset "v3, CS" begin
+        fmu = loadFMU(joinpath(FMU_DIR, "StateSpace.fmu"); type = :CS)
+        @named sspace = MTK.FMIComponent(
+            Val(3); fmu, communication_step_size = 1e-4, type = :CS,
+            reinitializealg = BrownFullBasicInit())
+        @test MTK.isinput(sspace.u)
+        @test MTK.isoutput(sspace.y)
+        @test !MTK.isinput(sspace.x) && !MTK.isoutput(sspace.x)
+
+        sys, prob = build_sspace_model(sspace)
+        sol = solve(prob, Rodas5P(autodiff = false); abstol = 1e-8, reltol = 1e-8)
+        @test SciMLBase.successful_retcode(sol)
+
+        @test truesol(
+            sol.t; idxs = [truesys.u, truesys.x, truesys.y]).u≈sol[[sys.u, sys.x, sys.y]] rtol=1e-2
+        @test truesol(sol.t; idxs = truesys.sspace.x).u≈sol(sol.t; idxs = sys.sspace.x).u rtol=1e-2
+    end
+end
+
+@testset "FMUs in a loop" begin
+    function build_looped_adders(adder1, adder2)
+        @variables x(t) = 1
+        @mtkbuild sys = ODESystem(
+            [D(x) ~ x, adder1.a ~ adder2.out2,
+                adder2.a ~ adder1.out2, adder1.b ~ 1.0, adder2.b ~ 2.0],
+            t;
+            systems = [adder1, adder2])
+        prob = ODEProblem(
+            sys, [adder1.c => 1.0, adder2.c => 1.0, adder1.a => 2.0], (0.0, 1.0))
+        return sys, prob
+    end
+    @named adder1 = SimpleAdder()
+    @named adder2 = SimpleAdder()
+    truesys, trueprob = build_looped_adders(adder1, adder2)
+    truesol = solve(trueprob, Tsit5(), reltol = 1e-8)
+    @test SciMLBase.successful_retcode(truesol)
+
+    @testset "v2, ME" begin
+        fmu = loadFMU(joinpath(FMU_DIR, "SimpleAdder.fmu"); type = :ME)
+        @named adder1 = MTK.FMIComponent(Val(2); fmu, type = :ME)
+        @named adder2 = MTK.FMIComponent(Val(2); fmu, type = :ME)
+        sys, prob = build_looped_adders(adder1, adder2)
+        sol = solve(prob, Rodas5P(autodiff = false); reltol = 1e-8)
+        @test SciMLBase.successful_retcode(sol)
+        @test truesol(sol.t;
+            idxs = [truesys.adder1.c, truesys.adder2.c]).u≈sol(
+            sol.t; idxs = [sys.adder1.c, sys.adder2.c]).u rtol=1e-7
+    end
+    @testset "v2, CS" begin
+        fmu = loadFMU(joinpath(FMU_DIR, "SimpleAdder.fmu"); type = :CS)
+        @named adder1 = MTK.FMIComponent(
+            Val(2); fmu, type = :CS, communication_step_size = 1e-3)
+        @named adder2 = MTK.FMIComponent(
+            Val(2); fmu, type = :CS, communication_step_size = 1e-3)
+        sys, prob = build_looped_adders(adder1, adder2)
+        sol = solve(prob, Tsit5(); reltol = 1e-8, initializealg = SciMLBase.OverrideInit(nlsolve = FastShortcutNLLSPolyalg(autodiff = AutoFiniteDiff())))
+        @test truesol(sol.t;
+            idxs = [truesys.adder1.c, truesys.adder2.c]).u≈sol(
+            sol.t; idxs = [sys.adder1.c, sys.adder2.c]).u rtol=1e-3
+    end
+
+    function build_looped_sspace(sspace1, sspace2)
+        @variables x(t) = 1
+        @mtkbuild sys = ODESystem([D(x) ~ x, sspace1.u ~ sspace2.x, sspace2.u ~ sspace1.y],
+            t; systems = [sspace1, sspace2])
+        prob = ODEProblem(sys, [sspace1.x => 1.0, sspace2.x => 1.0], (0.0, 1.0))
+        return sys, prob
+    end
+    @named sspace1 = StateSpace()
+    @named sspace2 = StateSpace()
+    truesys, trueprob = build_looped_sspace(sspace1, sspace2)
+    truesol = solve(trueprob, Rodas5P(), reltol = 1e-8)
+    @test SciMLBase.successful_retcode(truesol)
+
+    @testset "v3, ME" begin
+        fmu = loadFMU(joinpath(FMU_DIR, "StateSpace.fmu"); type = :ME)
+        @named sspace1 = MTK.FMIComponent(Val(3); fmu, type = :ME)
+        @named sspace2 = MTK.FMIComponent(Val(3); fmu, type = :ME)
+        sys, prob = build_looped_sspace(sspace1, sspace2)
+        sol = solve(prob, Rodas5P(autodiff = false); reltol = 1e-8)
+        @test SciMLBase.successful_retcode(sol)
+        @test truesol(sol.t;
+            idxs = [truesys.sspace1.x, truesys.sspace2.x]).u≈sol(
+            sol.t; idxs = [sys.sspace1.x, sys.sspace2.x]).u rtol=1e-7
+    end
+
+    @testset "v3, CS" begin
+        fmu = loadFMU(joinpath(FMU_DIR, "StateSpace.fmu"); type = :CS)
+        @named sspace1 = MTK.FMIComponent(
+            Val(3); fmu, type = :CS, communication_step_size = 1e-4)
+        @named sspace2 = MTK.FMIComponent(
+            Val(3); fmu, type = :CS, communication_step_size = 1e-4)
+        sys, prob = build_looped_sspace(sspace1, sspace2)
+        sol = solve(prob, Rodas5P(autodiff = false); reltol = 1e-8)
+        @test SciMLBase.successful_retcode(sol)
+        @test truesol(sol.t;
+            idxs = [truesys.sspace1.x, truesys.sspace2.x]).u≈sol(
+            sol.t; idxs = [sys.sspace1.x, sys.sspace2.x]).u rtol=1e-2
+    end
+end
diff --git a/test/fmi/fmus/SimpleAdder.fmu b/test/fmi/fmus/SimpleAdder.fmu
new file mode 100644
index 0000000000..ed12961bfb
Binary files /dev/null and b/test/fmi/fmus/SimpleAdder.fmu differ
diff --git a/test/fmi/fmus/SimpleAdder/README.md b/test/fmi/fmus/SimpleAdder/README.md
new file mode 100644
index 0000000000..b26c247bce
--- /dev/null
+++ b/test/fmi/fmus/SimpleAdder/README.md
@@ -0,0 +1,3 @@
+SimpleAdder.fmu is a v2 FMU built using OpenModelica. The file `SimpleAdder.mo` contains the
+modelica source code for the model. The file `buildFMU.mos` is the OpenModelica script used
+to build the FMU.
diff --git a/test/fmi/fmus/SimpleAdder/SimpleAdder.mo b/test/fmi/fmus/SimpleAdder/SimpleAdder.mo
new file mode 100644
index 0000000000..ce97f09a86
--- /dev/null
+++ b/test/fmi/fmus/SimpleAdder/SimpleAdder.mo
@@ -0,0 +1,12 @@
+block SimpleAdder
+	parameter Real value = 1.0;
+	input Real a;
+	input Real b;
+	Real c(start = 1.0, fixed = true);
+	output Real out;
+	output Real out2;
+equation
+	der(c) = out;
+	out = a + b + value;
+	out2 = 2 * c;
+end SimpleAdder;
diff --git a/test/fmi/fmus/SimpleAdder/buildFmu.mos b/test/fmi/fmus/SimpleAdder/buildFmu.mos
new file mode 100644
index 0000000000..df2c9bd58e
--- /dev/null
+++ b/test/fmi/fmus/SimpleAdder/buildFmu.mos
@@ -0,0 +1,12 @@
+OpenModelica.Scripting.loadFile("SimpleAdder.mo"); getErrorString();
+
+installPackage(Modelica);
+
+setCommandLineOptions("-d=newInst"); getErrorString();
+setCommandLineOptions("-d=initialization"); getErrorString();
+setCommandLineOptions("-d=-disableDirectionalDerivatives"); getErrorString();
+
+cd("output"); getErrorString();
+buildModelFMU(SimpleAdder, version = "2.0", fmuType = "me_cs"); getErrorString();
+system("unzip -l SimpleAdder.fmu | egrep -v 'sources|files' | tail -n+3 | grep -o '[A-Za-z._0-9/]*$' > BB.log")
+
diff --git a/test/fmi/fmus/StateSpace.fmu b/test/fmi/fmus/StateSpace.fmu
new file mode 100644
index 0000000000..6138e1fdb2
Binary files /dev/null and b/test/fmi/fmus/StateSpace.fmu differ
diff --git a/test/fmi/fmus/StateSpace/0001-tmp-commit.patch b/test/fmi/fmus/StateSpace/0001-tmp-commit.patch
new file mode 100644
index 0000000000..951af2d4fd
--- /dev/null
+++ b/test/fmi/fmus/StateSpace/0001-tmp-commit.patch
@@ -0,0 +1,425 @@
+From 10ac73996a7c09772ff31ccd6e030112d3bb5c43 Mon Sep 17 00:00:00 2001
+From: Aayush Sabharwal <aayush.sabharwal@gmail.com>
+Date: Wed, 8 Jan 2025 11:04:01 +0000
+Subject: tmp commit
+
+---
+ StateSpace/FMI3.xml |  31 ++----
+ StateSpace/config.h |  18 ++--
+ StateSpace/model.c  | 230 +++++++++++++-------------------------------
+ 3 files changed, 83 insertions(+), 196 deletions(-)
+
+diff --git a/StateSpace/FMI3.xml b/StateSpace/FMI3.xml
+index 8d2e4d9..002fbab 100644
+--- a/StateSpace/FMI3.xml
++++ b/StateSpace/FMI3.xml
+@@ -33,36 +33,23 @@
+         <UInt64 name="m" valueReference="1" description="Number of inputs" causality="structuralParameter" variability="tunable" start="3" min="0" max="5"/>
+         <UInt64 name="n" valueReference="2" description="Number of states" causality="structuralParameter" variability="tunable" start="3" min="0" max="5"/>
+         <UInt64 name="r" valueReference="3" description="Number of outputs" causality="structuralParameter" variability="tunable" start="3" min="0" max="5"/>
+-        <Float64 name="A" valueReference="4" description="Matrix coefficient A" causality="parameter" variability="tunable" start="1 0 0 0 1 0 0 0 1">
+-            <Dimension valueReference="2"/>
+-            <Dimension valueReference="2"/>
++        <Float64 name="A" valueReference="4" description="coefficient A" causality="parameter" variability="tunable" start="1">
+         </Float64>
+-        <Float64 name="B" valueReference="5" description="Matrix coefficient B" causality="parameter" variability="tunable" start="1 0 0 0 1 0 0 0 1">
+-            <Dimension valueReference="2"/>
+-            <Dimension valueReference="1"/>
++        <Float64 name="B" valueReference="5" description="coefficient B" causality="parameter" variability="tunable" start="1">
+         </Float64>
+-        <Float64 name="C" valueReference="6" description="Matrix coefficient C" causality="parameter" variability="tunable" start="1 0 0 0 1 0 0 0 1">
+-            <Dimension valueReference="3"/>
+-            <Dimension valueReference="2"/>
++        <Float64 name="C" valueReference="6" description="coefficient C" causality="parameter" variability="tunable" start="1">
+         </Float64>
+-        <Float64 name="D" valueReference="7" description="Matrix coefficient D" causality="parameter" variability="tunable" start="1 0 0 0 1 0 0 0 1">
+-            <Dimension valueReference="3"/>
+-            <Dimension valueReference="1"/>
++        <Float64 name="D" valueReference="7" description="coefficient D" causality="parameter" variability="tunable" start="1">
+         </Float64>
+-        <Float64 name="x0" valueReference="8" description="Initial state vector" causality="parameter" variability="tunable" start="0 0 0">
+-            <Dimension valueReference="2"/>
++        <Float64 name="x0" valueReference="8" description="Initial state" causality="parameter" variability="tunable" start="0">
+         </Float64>
+-        <Float64 name="u" valueReference="9" description="Input vector" causality="input" start="1 2 3">
+-            <Dimension valueReference="1"/>
++        <Float64 name="u" valueReference="9" description="Input" causality="input" start="1">
+         </Float64>
+-        <Float64 name="y" valueReference="10" description="Output vector" causality="output">
+-            <Dimension valueReference="3"/>
++        <Float64 name="y" valueReference="10" description="Output" causality="output">
+         </Float64>
+-        <Float64 name="x" valueReference="11" description="State vector" causality="local">
+-            <Dimension valueReference="2"/>
++        <Float64 name="x" valueReference="11" description="State" causality="local" initial = "exact" start = "1">
+         </Float64>
+-        <Float64 name="der(x)" valueReference="12" description="Derivative of the state vector" causality="local" derivative="11">
+-            <Dimension valueReference="2"/>
++        <Float64 name="der(x)" valueReference="12" description="Derivative of the state" causality="local" derivative="11">
+         </Float64>
+     </ModelVariables>
+ 
+diff --git a/StateSpace/config.h b/StateSpace/config.h
+index 707e500..8b422e2 100644
+--- a/StateSpace/config.h
++++ b/StateSpace/config.h
+@@ -44,15 +44,15 @@ typedef struct {
+     uint64_t m;
+     uint64_t n;
+     uint64_t r;
+-    double A[M_MAX][N_MAX];
+-    double B[M_MAX][N_MAX];
+-    double C[M_MAX][N_MAX];
+-    double D[M_MAX][N_MAX];
+-    double x0[N_MAX];
+-    double u[N_MAX];
+-    double y[N_MAX];
+-    double x[N_MAX];
+-    double der_x[N_MAX];
++    double A;
++    double B;
++    double C;
++    double D;
++    double x0;
++    double u;
++    double y;
++    double x;
++    double der_x;
+ } ModelData;
+ 
+ #endif /* config_h */
+diff --git a/StateSpace/model.c b/StateSpace/model.c
+index 8a47e74..9c170d8 100644
+--- a/StateSpace/model.c
++++ b/StateSpace/model.c
+@@ -3,77 +3,29 @@
+ 
+ 
+ void setStartValues(ModelInstance *comp) {
+-
+-    M(m) = 3;
+-    M(n) = 3;
+-    M(r) = 3;
+-
+     // identity matrix
+-    for (int i = 0; i < M_MAX; i++)
+-    for (int j = 0; j < N_MAX; j++) {
+-        M(A)[i][j] = i == j ? 1 : 0;
+-        M(B)[i][j] = i == j ? 1 : 0;
+-        M(C)[i][j] = i == j ? 1 : 0;
+-        M(D)[i][j] = i == j ? 1 : 0;
+-    }
+-
+-    for (int i = 0; i < M_MAX; i++) {
+-        M(u)[i] = i + 1;
+-    }
+-
+-    for (int i = 0; i < N_MAX; i++) {
+-        M(y)[i] = 0;
+-    }
+-
+-    for (int i = 0; i < N_MAX; i++) {
+-        M(x)[i] = M(x0)[i];
+-        M(x)[i] = 0;
+-    }
+-
++    M(m) = 1;
++    M(n) = 1;
++    M(r) = 1;
++
++    M(A) = 1;
++    M(B) = 1;
++    M(C) = 1;
++    M(D) = 1;
++
++    M(u) = 1;
++    M(y) = 0;
++    M(x) = 0;
++    M(x0) = 0;
+ }
+ 
+ Status calculateValues(ModelInstance *comp) {
+-
+-    // der(x) = Ax + Bu
+-    for (size_t i = 0; i < M(n); i++) {
+-
+-        M(der_x)[i] = 0;
+-
+-        for (size_t j = 0; j < M(n); j++) {
+-            M(der_x)[i] += M(A)[i][j] * M(x)[j];
+-        }
+-    }
+-
+-    for (size_t i = 0; i < M(n); i++) {
+-
+-        for (size_t j = 0; j < M(r); j++) {
+-            M(der_x)[i] += M(B)[i][j] * M(u)[j];
+-        }
+-    }
+-
+-
+-    // y = Cx + Du
+-    for (size_t i = 0; i < M(r); i++) {
+-
+-        M(y)[i] = 0;
+-
+-        for (size_t j = 0; j < M(n); j++) {
+-            M(y)[i] += M(C)[i][j] * M(x)[j];
+-        }
+-    }
+-
+-    for (size_t i = 0; i < M(r); i++) {
+-
+-        for (size_t j = 0; j < M(m); j++) {
+-            M(y)[i] += M(D)[i][j] * M(u)[j];
+-        }
+-    }
+-
++    M(der_x) = M(A) * M(x) + M(B) * M(u);
++    M(y) = M(C) * M(x) + M(D) * M(u);
+     return OK;
+ }
+ 
+ Status getFloat64(ModelInstance* comp, ValueReference vr, double values[], size_t nValues, size_t* index) {
+-
+     calculateValues(comp);
+ 
+     switch (vr) {
+@@ -82,66 +34,40 @@ Status getFloat64(ModelInstance* comp, ValueReference vr, double values[], size_
+             values[(*index)++] = comp->time;
+             return OK;
+         case vr_A:
+-            ASSERT_NVALUES((size_t)(M(n) * M(n)));
+-            for (size_t i = 0; i < M(n); i++) {
+-                for (size_t j = 0; j < M(n); j++) {
+-                    values[(*index)++] = M(A)[i][j];
+-                }
+-            }
++            ASSERT_NVALUES(1);
++            values[(*index)++] = M(A);
+             return OK;
+         case vr_B:
+-            ASSERT_NVALUES((size_t)(M(m) * M(n)));
+-            for (size_t i = 0; i < M(m); i++) {
+-                for (size_t j = 0; j < M(n); j++) {
+-                    values[(*index)++] = M(B)[i][j];
+-                }
+-            }
++            ASSERT_NVALUES(1);
++            values[(*index)++] = M(B);
+             return OK;
+         case vr_C:
+-            ASSERT_NVALUES((size_t)(M(r) * M(n)));
+-            for (size_t i = 0; i < M(r); i++) {
+-                for (size_t j = 0; j < M(n); j++) {
+-                    values[(*index)++] = M(C)[i][j];
+-                }
+-            }
++            ASSERT_NVALUES(1);
++            values[(*index)++] = M(C);
+             return OK;
+         case vr_D:
+-            ASSERT_NVALUES((size_t)(M(r) * M(m)));
+-            for (size_t i = 0; i < M(r); i++) {
+-                for (size_t j = 0; j < M(m); j++) {
+-                    values[(*index)++] = M(D)[i][j];
+-                }
+-            }
++            ASSERT_NVALUES(1);
++            values[(*index)++] = M(D);
+             return OK;
+         case vr_x0:
+-            ASSERT_NVALUES((size_t)M(n));
+-            for (size_t i = 0; i < M(n); i++) {
+-                values[(*index)++] = M(x0)[i];
+-            }
++            ASSERT_NVALUES(1);
++            values[(*index)++] = M(x0);
+             return OK;
+         case vr_u:
+-            ASSERT_NVALUES((size_t)M(m));
+-            for (size_t i = 0; i < M(m); i++) {
+-                values[(*index)++] = M(u)[i];
+-            }
++            ASSERT_NVALUES(1);
++            values[(*index)++] = M(u);
+             return OK;
+         case vr_y:
+-            ASSERT_NVALUES((size_t)M(r));
+-            for (size_t i = 0; i < M(r); i++) {
+-                values[(*index)++] = M(y)[i];
+-            }
++            ASSERT_NVALUES(1);
++            values[(*index)++] = M(y);
+             return OK;
+         case vr_x:
+-            ASSERT_NVALUES((size_t)M(n));
+-            for (size_t i = 0; i < M(n); i++) {
+-                values[(*index)++] = M(x)[i];
+-            }
++            ASSERT_NVALUES(1);
++            values[(*index)++] = M(x);
+             return OK;
+         case vr_der_x:
+-            ASSERT_NVALUES((size_t)M(n));
+-            for (size_t i = 0; i < M(n); i++) {
+-                values[(*index)++] = M(der_x)[i];
+-            }
++            ASSERT_NVALUES(1);
++            values[(*index)++] = M(der_x);
+             return OK;
+         default:
+             logError(comp, "Get Float64 is not allowed for value reference %u.", vr);
+@@ -153,58 +79,40 @@ Status setFloat64(ModelInstance* comp, ValueReference vr, const double values[],
+ 
+     switch (vr) {
+     case vr_A:
+-        ASSERT_NVALUES((size_t)(M(n) * M(n)));
+-        for (size_t i = 0; i < M(n); i++) {
+-            for (size_t j = 0; j < M(n); j++) {
+-                M(A)[i][j] = values[(*index)++];
+-            }
+-        }
++        ASSERT_NVALUES(1);
++        M(A) = values[(*index)++];
+         break;
+     case vr_B:
+-        ASSERT_NVALUES((size_t)(M(n) * M(m)));
+-        for (size_t i = 0; i < M(n); i++) {
+-            for (size_t j = 0; j < M(m); j++) {
+-                M(B)[i][j] = values[(*index)++];
+-            }
+-        }
++        ASSERT_NVALUES(1);
++        M(B) = values[(*index)++];
+         break;
+     case vr_C:
+-        ASSERT_NVALUES((size_t)(M(r) * M(n)));
+-        for (size_t i = 0; i < M(r); i++) {
+-            for (size_t j = 0; j < M(n); j++) {
+-                M(C)[i][j] = values[(*index)++];
+-            }
+-        }
++        ASSERT_NVALUES(1);
++        M(C) = values[(*index)++];
+         break;
+     case vr_D:
+-        ASSERT_NVALUES((size_t)(M(r) * M(m)));
+-        for (size_t i = 0; i < M(r); i++) {
+-            for (size_t j = 0; j < M(m); j++) {
+-                M(D)[i][j] = values[(*index)++];
+-            }
+-        }
++        ASSERT_NVALUES(1);
++        M(D) = values[(*index)++];
+         break;
+     case vr_x0:
+-        ASSERT_NVALUES((size_t)M(n));
+-        for (size_t i = 0; i < M(n); i++) {
+-            M(x0)[i] = values[(*index)++];
+-        }
++        ASSERT_NVALUES(1);
++        M(x0) = values[(*index)++];
+         break;
+     case vr_u:
+-        ASSERT_NVALUES((size_t)M(m));
+-        for (size_t i = 0; i < M(m); i++) {
+-            M(u)[i] = values[(*index)++];
+-        }
++        ASSERT_NVALUES(1);
++        M(u) = values[(*index)++];
++        break;
++    case vr_y:
++        ASSERT_NVALUES(1);
++        M(y) = values[(*index)++];
+         break;
+     case vr_x:
+-        if (comp->state != ContinuousTimeMode && comp->state != EventMode) {
+-            logError(comp, "Variable \"x\" can only be set in Continuous Time Mode and Event Mode.");
+-            return Error;
+-        }
+-        ASSERT_NVALUES((size_t)M(n));
+-        for (size_t i = 0; i < M(n); i++) {
+-            M(x)[i] = values[(*index)++];
+-        }
++        ASSERT_NVALUES(1);
++        M(x) = values[(*index)++];
++        break;
++    case vr_der_x:
++        ASSERT_NVALUES(1);
++        M(der_x) = values[(*index)++];
+         break;
+     default:
+         logError(comp, "Set Float64 is not allowed for value reference %u.", vr);
+@@ -217,8 +125,6 @@ Status setFloat64(ModelInstance* comp, ValueReference vr, const double values[],
+ }
+ 
+ Status getUInt64(ModelInstance* comp, ValueReference vr, uint64_t values[], size_t nValues, size_t* index) {
+-
+-
+     switch (vr) {
+         case vr_m:
+             ASSERT_NVALUES(1);
+@@ -289,49 +195,43 @@ Status eventUpdate(ModelInstance *comp) {
+ }
+ 
+ size_t getNumberOfContinuousStates(ModelInstance* comp) {
+-    return (size_t)M(n);
++    return M(n) == 1 ? M(n) : 1;
+ }
+ 
+ Status getContinuousStates(ModelInstance* comp, double x[], size_t nx) {
+ 
+-    if (nx != M(n)) {
+-        logError(comp, "Expected nx=%zu but was %zu.", M(n), nx);
++    if (nx != 1) {
++        logError(comp, "Expected nx=%zu but was %zu.", 1, nx);
+         return Error;
+     }
+ 
+-    for (size_t i = 0; i < M(n); i++) {
+-        x[i] = M(x)[i];
+-    }
++    x[0] = M(x);
+ 
+     return OK;
+ }
+ 
+ Status setContinuousStates(ModelInstance* comp, const double x[], size_t nx) {
+ 
+-    if (nx != M(n)) {
+-        logError(comp, "Expected nx=%zu but was %zu.", M(n), nx);
++    if (nx != 1) {
++        logError(comp, "Expected nx=%zu but was %zu.", 1, nx);
+         return Error;
+     }
+ 
+-    for (size_t i = 0; i < M(n); i++) {
+-        M(x)[i] = x[i];
+-    }
++    M(x) = x[0];
+ 
+     return OK;
+ }
+ 
+ Status getDerivatives(ModelInstance* comp, double dx[], size_t nx) {
+ 
+-    if (nx != M(n)) {
+-        logError(comp, "Expected nx=%zu but was %zu.", M(n), nx);
++    if (nx != 1) {
++        logError(comp, "Expected nx=%zu but was %zu.", 1, nx);
+         return Error;
+     }
+ 
+     calculateValues(comp);
+ 
+-    for (size_t i = 0; i < M(n); i++) {
+-        dx[i] = M(der_x)[i];
+-    }
++    dx[0] = M(der_x);
+ 
+     return OK;
+ }
+-- 
+2.34.1
+
diff --git a/test/fmi/fmus/StateSpace/README.md b/test/fmi/fmus/StateSpace/README.md
new file mode 100644
index 0000000000..108147443d
--- /dev/null
+++ b/test/fmi/fmus/StateSpace/README.md
@@ -0,0 +1,7 @@
+StateSpace.fmu is a v3 FMU built using a modified version of the StateSpace FMU in Modelica's
+Reference-FMUs. The link to the specific commit hash used is
+
+[https://github.com/modelica/Reference-FMUs/tree/0e1374cad0a7f0583a583ce189060ba7a67c27a7]()
+
+The git patch containing the changes is provided as `0001-tmp-commit.patch`. The FMU can be
+built using the instructions provided in the repository's README.
diff --git a/test/initial_values.jl b/test/initial_values.jl
index a2931e22b6..c31eca33f8 100644
--- a/test/initial_values.jl
+++ b/test/initial_values.jl
@@ -174,3 +174,11 @@ end
     @test_throws ModelingToolkit.UnexpectedSymbolicValueInVarmap ODEProblem(
         sys, [x => 1, y => 2], (0.0, 1.0), [p => 2q, q => 3p])
 end
+
+@testset "`add_fallbacks!` checks scalarized array parameters correctly" begin
+    @variables x(t)[1:2]
+    @parameters p[1:2, 1:2]
+    @mtkbuild sys = ODESystem(D(x) ~ p * x, t)
+    # used to throw a `MethodError` complaining about `getindex(::Nothing, ::CartesianIndex{2})`
+    @test_throws ModelingToolkit.MissingParametersError ODEProblem(sys, [x => ones(2)], (0.0, 1.0))
+end
diff --git a/test/initializationsystem.jl b/test/initializationsystem.jl
index 0a4ef18038..def45a7abf 100644
--- a/test/initializationsystem.jl
+++ b/test/initializationsystem.jl
@@ -792,10 +792,10 @@ end
         @variables x=1.0 y=3.0
         @parameters p q
 
-        @mtkbuild sys = NonlinearSystem(
-            [(x - p)^2 + (y - q)^3 ~ 0, x - q ~ 0]; defaults = [q => missing],
+        @named sys = NonlinearSystem(
+            [x + y - p ~ 0, x - q ~ 0]; defaults = [q => missing],
             guesses = [q => 1.0], initialization_eqs = [p^2 + q^2 + 2p * q ~ 0])
-
+        sys = complete(sys)
         for (probT, algs) in prob_alg_combinations
             prob = probT(sys, [], [p => 2.0])
             @test prob.f.initialization_data !== nothing
diff --git a/test/odesystem.jl b/test/odesystem.jl
index c1b0372c0a..a635c3dad9 100644
--- a/test/odesystem.jl
+++ b/test/odesystem.jl
@@ -1544,6 +1544,32 @@ end
     @test all(x -> any(isapprox(x, atol = 1e-6), sol2.t), expected_tstops)
 end
 
+# Test `isequal`
+@testset "`isequal`" begin
+    @variables X(t)
+    @parameters p d
+    eq = D(X) ~ p - d*X
+
+    osys1 = complete(ODESystem([eq], t; name = :osys))
+    osys2 = complete(ODESystem([eq], t; name = :osys))
+    @test osys1 == osys2 # true
+
+    continuous_events = [[X ~ 1.0] => [X ~ X + 5.0]]
+    discrete_events = [5.0 => [d ~ d / 2.0]]
+
+    osys1 = complete(ODESystem([eq], t; name = :osys, continuous_events))
+    osys2 = complete(ODESystem([eq], t; name = :osys))
+    @test osys1 !== osys2 
+
+    osys1 = complete(ODESystem([eq], t; name = :osys, discrete_events))
+    osys2 = complete(ODESystem([eq], t; name = :osys))
+    @test osys1 !== osys2
+
+    osys1 = complete(ODESystem([eq], t; name = :osys, continuous_events))
+    osys2 = complete(ODESystem([eq], t; name = :osys, discrete_events))
+    @test osys1 !== osys2 
+end
+
 @testset "dae_order_lowering basic test" begin
     @parameters a
     @variables x(t) y(t) z(t)
@@ -1599,4 +1625,4 @@ end
 
     prob = ODEProblem{false}(lowered_dae_sys; u0_constructor = x -> SVector(x...))
     @test prob.u0 isa SVector
-end
\ No newline at end of file
+end
diff --git a/test/runtests.jl b/test/runtests.jl
index 95154e550e..52875fdae5 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -4,6 +4,12 @@ import REPL
 
 const GROUP = get(ENV, "GROUP", "All")
 
+function activate_fmi_env()
+    Pkg.activate("fmi")
+    Pkg.develop(PackageSpec(path = dirname(@__DIR__)))
+    Pkg.instantiate()
+end
+
 function activate_extensions_env()
     Pkg.activate("extensions")
     Pkg.develop(PackageSpec(path = dirname(@__DIR__)))
@@ -115,6 +121,11 @@ end
         @safetestset "Analysis Points Test" include("downstream/analysis_points.jl")
     end
 
+    if GROUP == "All" || GROUP == "FMI"
+        activate_fmi_env()
+        @safetestset "FMI Extension Test" include("fmi/fmi.jl")
+    end
+
     if GROUP == "All" || GROUP == "Extensions"
         activate_extensions_env()
         @safetestset "HomotopyContinuation Extension Test" include("extensions/homotopy_continuation.jl")
diff --git a/test/sdesystem.jl b/test/sdesystem.jl
index 94f3e9cbfc..ade91b60b9 100644
--- a/test/sdesystem.jl
+++ b/test/sdesystem.jl
@@ -869,6 +869,73 @@ end
     @test length(observed(sys)) == 1
 end
 
+@testset "SDEFunctionExpr" begin
+    @parameters σ ρ β
+    @variables x(tt) y(tt) z(tt)
+
+    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]
+
+    @named sys = ODESystem(eqs, tt, [x, y, z], [σ, ρ, β])
+    
+    @named de = SDESystem(eqs, noiseeqs, tt, [x, y, z], [σ, ρ, β], tspan = (0.0, 10.0))
+    de = complete(de)
+    
+    f = SDEFunctionExpr(de)
+    @test f isa Expr
+
+    @testset "Configuration Tests" begin
+        # Test with `tgrad`
+        f_tgrad = SDEFunctionExpr(de; tgrad = true)
+        @test f_tgrad isa Expr
+
+        # Test with `jac`
+        f_jac = SDEFunctionExpr(de; jac = true)
+        @test f_jac isa Expr
+
+        # Test with sparse Jacobian
+        f_sparse = SDEFunctionExpr(de; sparse = true)
+        @test f_sparse isa Expr
+    end
+
+    @testset "Ordering Tests" begin
+        dvs = [z, y, x]
+        ps = [β, ρ, σ]
+        f_order = SDEFunctionExpr(de, dvs, ps)
+        @test f_order isa Expr
+    end
+end
+
+@testset "SDESystem Equality with events" begin
+    @variables X(t)
+    @parameters p d
+    @brownian a
+    seq = D(X) ~ p - d*X + a
+    @mtkbuild ssys1 = System([seq], t; name = :ssys)
+    @mtkbuild ssys2 = System([seq], t; name = :ssys)
+    @test ssys1 == ssys2 # true
+
+    continuous_events = [[X ~ 1.0] => [X ~ X + 5.0]]
+    discrete_events = [5.0 => [d ~ d / 2.0]]
+
+    @mtkbuild ssys1 = System([seq], t; name = :ssys, continuous_events)
+    @mtkbuild ssys2 = System([seq], t; name = :ssys)
+    @test ssys1 !== ssys2 
+
+    @mtkbuild ssys1 = System([seq], t; name = :ssys, discrete_events)
+    @mtkbuild ssys2 = System([seq], t; name = :ssys)
+    @test ssys1 !== ssys2
+
+    @mtkbuild ssys1 = System([seq], t; name = :ssys, continuous_events)
+    @mtkbuild ssys2 = System([seq], t; name = :ssys, discrete_events)
+    @test ssys1 !== ssys2 
+end
+
 @testset "Error when constructing SDESystem without `structural_simplify`" begin
     @parameters σ ρ β
     @variables x(tt) y(tt) z(tt)
@@ -886,4 +953,4 @@ end
     @test_throws ErrorException("SDESystem constructed by defining Brownian variables with @brownian must be simplified by calling `structural_simplify` before a SDEProblem can be constructed.") SDEProblem(de, u0map, (0.0, 100.0), parammap)
     de = structural_simplify(de)
     @test SDEProblem(de, u0map, (0.0, 100.0), parammap) isa SDEProblem
-end
+end
\ No newline at end of file
diff --git a/test/structural_transformation/utils.jl b/test/structural_transformation/utils.jl
index 863e091aad..621aa8b8a8 100644
--- a/test/structural_transformation/utils.jl
+++ b/test/structural_transformation/utils.jl
@@ -92,7 +92,7 @@ end
         @mtkbuild sys = ODESystem(
             [D(y) ~ foo(x), D(x) ~ sum(y), zeros(2) ~ foo(prod(z))], t)
         @test length(equations(sys)) == 5
-        @test length(observed(sys)) == 2
+        @test length(observed(sys)) == 4
         prob = ODEProblem(
             sys, [y => ones(2), z => 2ones(2), x => 3.0], (0.0, 1.0), [foo => _tmp_fn2])
         @test_nowarn prob.f(prob.u0, prob.p, 0.0)