Error trying to chain 2 simple ODEs in ModelingToolkit

simpleconnectiontest.jl

using ModelingToolkit

@parameters t a
@variables x(t) f(t)
D = Differential(t)

eqs = [
  D(x) ~ -a*x + f
]
decay = name -> ODESystem(eqs, name=name)

decay1 = decay(:decay1)
decay2 = decay(:decay2)

connected = ODESystem([
  decay2.f ~ decay1.x
], t, systems=[decay1, decay2])

reduced_system = alias_elimination(connected; conservative=false)
@show reduced_system

x0 = [
  decay1.x => 1.0
  decay1.f => 0.0
  decay2.x => 1.0
]
p = [
  decay1.a => 0.1
  decay2.a => 0.2
]

# produces ERROR: BoundsError: attempt to access 2×2 Matrix{Float64} at index [2, 3]
prob = ODEProblem(reduced_system, x0, (0.0, 20.0), p)
sol = solve(prob)

Full output:

julia> using ModelingToolkit

julia> @parameters t a
(t, a)

julia> @variables x(t) f(t)
(x(t), f(t))

julia> D = Differential(t)
(::Differential) (generic function with 2 methods)

julia> eqs = [
         D(x) ~ -a*x + f
         # D(f) ~ 0
       ]
1-element Vector{Equation}:
 Differential(t)(x(t)) ~ f(t) - (a*(x(t)))

julia> decay = name -> ODESystem(eqs, name=name)
#1 (generic function with 1 method)

julia> decay1 = decay(:decay1)
Model decay1 with 1 equations
States (2):
  x(t)
  f(t)
Parameters (1):
  a

julia> decay2 = decay(:decay2)
Model decay2 with 1 equations
States (2):
  x(t)
  f(t)
Parameters (1):
  a

julia> connected = ODESystem([
         decay2.f ~ decay1.x
       ], t, systems=[decay1, decay2])
Model ##ODESystem#257 with 3 equations
States (4):
  decay2₊f(t)
  decay1₊x(t)
  decay1₊f(t)
  decay2₊x(t)
Parameters (2):
  decay1₊a
  decay2₊a

julia> reduced_system = alias_elimination(connected; conservative=false)
Model ##ODESystem#258 with 2 equations
States (3):
  decay1₊x(t)
  decay1₊f(t)
  decay2₊x(t)
Parameters (2):
  decay1₊a
  decay2₊a
Incidence matrix:
 ×  ⋅  ×  ⋅  ×
 ×  ×  ⋅  ×  ⋅

julia> @show reduced_system
reduced_system = Model ##ODESystem#258 with 2 equations
States (3):
  decay1₊x(t)
  decay1₊f(t)
  decay2₊x(t)
Parameters (2):
  decay1₊a
  decay2₊a
Incidence matrix:
 ×  ⋅  ×  ⋅  ×
 ×  ×  ⋅  ×  ⋅
Model ##ODESystem#258 with 2 equations
States (3):
  decay1₊x(t)
  decay1₊f(t)
  decay2₊x(t)
Parameters (2):
  decay1₊a
  decay2₊a
Incidence matrix:
 ×  ⋅  ×  ⋅  ×
 ×  ×  ⋅  ×  ⋅

julia> x0 = [
         decay1.x => 1.0
         decay1.f => 0.0
         decay2.x => 1.0
         # decay2.f => 0.0,
       ]
3-element Vector{Pair{Term{Real}, Float64}}:
 decay1₊x(t) => 1.0
 decay1₊f(t) => 0.0
 decay2₊x(t) => 1.0

julia> p = [
         decay1.a => 0.1
         decay2.a => 0.2
       ]
2-element Vector{Pair{Sym{ModelingToolkit.Parameter{Real}}, Float64}}:
 decay1₊a => 0.1
 decay2₊a => 0.2

julia> prob = ODEProblem(reduced_system, x0, (0.0, 20.0), p)
ERROR: BoundsError: attempt to access 2×2 Matrix{Float64} at index [2, 3]
Stacktrace:
  [1] setindex!(::Matrix{Float64}, ::Int64, ::Int64, ::Int64)
    @ Base ./array.jl:841
  [2] calculate_massmatrix(sys::ODESystem; simplify::Bool)
    @ ModelingToolkit ~/.julia/packages/ModelingToolkit/pr3qG/src/systems/diffeqs/abstractodesystem.jl:88
  [3] calculate_massmatrix
    @ ~/.julia/packages/ModelingToolkit/pr3qG/src/systems/diffeqs/abstractodesystem.jl:82 [inlined]
  [4] (ODEFunction{true, F, TMM, Ta, Tt, TJ, JVP, VJP, JP, SP, TW, TWt, TPJ, S, S2, O, TCV} where {F, TMM, Ta, Tt, TJ, JVP, VJP, JP, SP, TW, TWt, TPJ, S, S2, O, TCV})(sys::ODESystem, dvs::Vector{Term{Real}}, ps::Vector{Sym{ModelingToolkit.Parameter{Real}}}, u0::Vector{Float64}; version::Nothing, tgrad::Bool, jac::Bool, eval_expression::Bool, sparse::Bool, simplify::Bool, eval_module::Module, kwargs::Base.Iterators.Pairs{Symbol, Any, Tuple{Symbol, Symbol, Symbol}, NamedTuple{(:checkbounds, :linenumbers, :parallel), Tuple{Bool, Bool, ModelingToolkit.SerialForm}}})
    @ ModelingToolkit ~/.julia/packages/ModelingToolkit/pr3qG/src/systems/diffeqs/abstractodesystem.jl:156
  [5] process_DEProblem(constructor::Type, sys::ODESystem, u0map::Vector{Pair{Term{Real}, Float64}}, parammap::Vector{Pair{Sym{ModelingToolkit.Parameter{Real}}, Float64}}; version::Nothing, tgrad::Bool, jac::Bool, checkbounds::Bool, sparse::Bool, simplify::Bool, linenumbers::Bool, parallel::ModelingToolkit.SerialForm, eval_expression::Bool, kwargs::Base.Iterators.Pairs{Union{}, Union{}, Tuple{}, NamedTuple{(), Tuple{}}})
    @ ModelingToolkit ~/.julia/packages/ModelingToolkit/pr3qG/src/systems/diffeqs/abstractodesystem.jl:262
  [6] process_DEProblem
    @ ~/.julia/packages/ModelingToolkit/pr3qG/src/systems/diffeqs/abstractodesystem.jl:245 [inlined]
  [7] (ODEProblem{true, tType, isinplace, P, F, K, PT} where {tType, isinplace, P, F, K, PT})(sys::ODESystem, u0map::Vector{Pair{Term{Real}, Float64}}, tspan::Tuple{Float64, Float64}, parammap::Vector{Pair{Sym{ModelingToolkit.Parameter{Real}}, Float64}}; kwargs::Base.Iterators.Pairs{Union{}, Union{}, Tuple{}, NamedTuple{(), Tuple{}}})
    @ ModelingToolkit ~/.julia/packages/ModelingToolkit/pr3qG/src/systems/diffeqs/abstractodesystem.jl:294
  [8] (ODEProblem{true, tType, isinplace, P, F, K, PT} where {tType, isinplace, P, F, K, PT})(sys::ODESystem, u0map::Vector{Pair{Term{Real}, Float64}}, tspan::Tuple{Float64, Float64}, parammap::Vector{Pair{Sym{ModelingToolkit.Parameter{Real}}, Float64}})
    @ ModelingToolkit ~/.julia/packages/ModelingToolkit/pr3qG/src/systems/diffeqs/abstractodesystem.jl:294
  [9] ODEProblem(::ODESystem, ::Vector{Pair{Term{Real}, Float64}}, ::Vararg{Any, N} where N; kwargs::Base.Iterators.Pairs{Union{}, Union{}, Tuple{}, NamedTuple{(), Tuple{}}})
    @ ModelingToolkit ~/.julia/packages/ModelingToolkit/pr3qG/src/systems/diffeqs/abstractodesystem.jl:274
 [10] ODEProblem(::ODESystem, ::Vector{Pair{Term{Real}, Float64}}, ::Vararg{Any, N} where N)
    @ ModelingToolkit ~/.julia/packages/ModelingToolkit/pr3qG/src/systems/diffeqs/abstractodesystem.jl:274
 [11] top-level scope
    @ REPL[14]:1