Forward MTK system to `ODEFunction` constructor in `ODEFunctionExpr`

[devmotion]

Is your feature request related to a problem? Please describe.

I'm surprised that ODEFunctionExpr(sys) returns an expression that only sets a subset of properties of sys, and in particular does not forward the system sys.

For instance, the example from the README gives

julia> using ModelingToolkit

julia> using ModelingToolkit: t_nounits as t, D_nounits as D

julia> @parameters σ ρ β
3-element Vector{Num}:
 σ
 ρ
 β

julia> @variables x(t) y(t) z(t)
3-element Vector{Num}:
 x(t)
 y(t)
 z(t)

julia> eqs = [D(D(x)) ~ σ * (y - x),
           D(y) ~ x * (ρ - z) - y,
           D(z) ~ x * y - β * z]
3-element Vector{Equation}:
 Differential(t)(Differential(t)(x(t))) ~ (-x(t) + y(t))*σ
 Differential(t)(y(t)) ~ -y(t) + x(t)*(-z(t) + ρ)
 Differential(t)(z(t)) ~ x(t)*y(t) - z(t)*β

julia> @mtkbuild sys = ODESystem(eqs, t)
Model sys with 4 equations
Unknowns (4):
  xˍt(t)
  y(t)
  z(t)
  x(t)
Parameters (3):
  σ
  ρ
  β
Incidence matrix:4×8 SparseArrays.SparseMatrixCSC{Num, Int64} with 13 stored entries:
 ⋅  ×  ⋅  ×  ×  ⋅  ⋅  ⋅
 ⋅  ×  ×  ×  ⋅  ⋅  ×  ⋅
 ⋅  ×  ×  ×  ⋅  ⋅  ⋅  ×
 ×  ⋅  ⋅  ⋅  ⋅  ×  ⋅  ⋅

julia> ex = ODEFunctionExpr(sys)
quote
    var"##f#226" = (ModelingToolkit.ODEFunctionClosure)(function (ˍ₋arg1, ˍ₋arg2, t)
                begin
                    begin
                        begin
                            (SymbolicUtils.Code.create_array)(typeof(ˍ₋arg1), nothing, Val{1}(), Val{(4,)}(), (*)((+)((*)(-1, ˍ₋arg1[4]), ˍ₋arg1[2]), ˍ₋arg2[3]), (+)((*)(-1, ˍ₋arg1[2]), (*)(ˍ₋arg1[4], (+)((*)(-1, ˍ₋arg1[3]), ˍ₋arg2[1]))), (+)((*)(ˍ₋arg1[4], ˍ₋arg1[2]), (*)((*)(-1, ˍ₋arg1[3]), ˍ₋arg2[2])), ˍ₋arg1[1])
                        end
                    end
                end
            end, function (ˍ₋out, ˍ₋arg1, ˍ₋arg2, t)
                begin
                    begin
                        begin
                            #= /Users/david/.julia/packages/SymbolicUtils/c0xQb/src/code.jl:422 =# @inbounds begin
                                    ˍ₋out[1] = (*)((+)((*)(-1, ˍ₋arg1[4]), ˍ₋arg1[2]), ˍ₋arg2[3])
                                    ˍ₋out[2] = (+)((*)(-1, ˍ₋arg1[2]), (*)(ˍ₋arg1[4], (+)((*)(-1, ˍ₋arg1[3]), ˍ₋arg2[1])))
                                    ˍ₋out[3] = (+)((*)(ˍ₋arg1[4], ˍ₋arg1[2]), (*)((*)(-1, ˍ₋arg1[3]), ˍ₋arg2[2]))
                                    ˍ₋out[4] = ˍ₋arg1[1]
                                    nothing
                                end
                        end
                    end
                end
            end)
    var"##tgrad#227" = nothing
    var"##jac#228" = nothing
    M = LinearAlgebra.UniformScaling{Bool}(true)
    ODEFunction{true}(var"##f#226", jac = var"##jac#228", tgrad = var"##tgrad#227", mass_matrix = M, jac_prototype = nothing, sparsity = nothing, observed = nothing)
end

Describe the solution you’d like

I'd like ODEFunctionExpr to return an ODEFunction expression that includes sys - similar to the definition of ODEFunction in https://github.com/SciML/ModelingToolkit.jl/blob/539fe81a800a8fc559234ba3cfa50720ab5256f1/src/systems/diffeqs/abstractodesystem.jl#L499. Maybe sys could be wrapped as an addition field :sys in the ODEFunctionClosure? Then one would not have to define an additional variable in the expression and the constructor in SciMLBase would still use it automatically: https://github.com/SciML/SciMLBase.jl/blob/4fdae7eb4a6f9bf22236bf2aabd0c210459b2503/src/scimlfunctions.jl#L2248

[ChrisRackauckas]

Yes the Expr form is missing the pass through of a few things and it's just not implemented yet. It wouldn't be difficult to add it though.