# Using ModelingToolkit.jl to model a matrix ODE

**URL:** <https://discourse.julialang.org/t/using-modelingtoolkit-jl-to-model-a-matrix-ode/82490>\
**Category:** Modelling & Simulations\
**Tags:** question, diffeq, modelingtoolkit\
**Created:** [June 9, 2022, 11:23am UTC](https://discourse.julialang.org/t/using-modelingtoolkit-jl-to-model-a-matrix-ode/82490 "2022-06-09T11:23:56Z")\
**Posts on this page:** 4\
**Page:** 1

<div class="post-metadata">

**Author:** ![golanor](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/golanor/32/108140_2.png) [@golanor](https://discourse.julialang.org/u/golanor)\
**Post date:** [June 9, 2022, 11:23am UTC](https://discourse.julialang.org/t/using-modelingtoolkit-jl-to-model-a-matrix-ode/82490/1 "2022-06-09T11:23:56Z")

</div>

Hi,

I’m trying to create a matrix ODE using ModelingToolkit,jl (to simplify autodiff and code optimization), but I’m receiving an error which I did not find in any of the issues of the library. Any idea what I’m doing wrong?

MWE:

```julia
using ModelingToolkit
using DifferentialEquations: solve

@variables x y[1:2, 1:2](x)
@parameters f g

mat = [1.0 2.0; 0.0 0.5]

D = Differential(x)
@named sys = ODESystem([D(y) ~ g * y + f * mat * y])

prob = ODEProblem(sys, [y => [1.0 0.0; 0.0 0.0]], [g => 1.0, f => 2.0])

solve(prob)

```

Stacktrace:

```julia
Stacktrace:
  [1] copyto!(dest::Vector{Any}, src::SymbolicUtils.Mul{Real, Int64, Dict{Any, Number}, Nothing})
    @ Base ./abstractarray.jl:890
  [2] _collect(cont::UnitRange{Int64}, itr::SymbolicUtils.Mul{Real, Int64, Dict{Any, Number}, Nothing}, #unused#::Base.HasEltype, isz::Base.HasLength)
    @ Base ./array.jl:655
  [3] collect(itr::SymbolicUtils.Mul{Real, Int64, Dict{Any, Number}, Nothing})
    @ Base ./array.jl:649
  [4] broadcastable(x::SymbolicUtils.Mul{Real, Int64, Dict{Any, Number}, Nothing})
    @ Base.Broadcast ./broadcast.jl:704
  [5] broadcasted
    @ ./broadcast.jl:1300 [inlined]
  [6] broadcast(::typeof(*), ::SymbolicUtils.Mul{Real, Int64, Dict{Any, Number}, Nothing}, ::Symbolics.ArrayOp{AbstractMatrix{Real}})
    @ Base.Broadcast ./broadcast.jl:798
  [7] materialize(bc::Base.Broadcast.Broadcasted{Symbolics.SymWrapBroadcast, Nothing, typeof(*), Tuple{Num, Symbolics.Arr{Num, 2}}})
    @ Symbolics ~/.julia/packages/Symbolics/H8dtg/src/array-lib.jl:150
  [8] broadcast_preserving_zero_d
    @ ./broadcast.jl:849 [inlined]
  [9] *(A::Num, B::Symbolics.Arr{Num, 2})
    @ Base ./arraymath.jl:52
 [10] *(::Num, ::Num, ::Symbolics.Arr{Num, 2})
    @ Base ./operators.jl:655
 [11] top-level scope
    @ ~/Code/julia/control/src/mwe.jl:10
 [12] include(fname::String)
    @ Base.MainInclude ./client.jl:451
 [13] top-level scope
    @ REPL[3]:1
in expression starting at ./mwe.jl:10

```

---

<div class="post-metadata">

**Author:** ![FlorianD](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/floriand/32/38696_2.png) [@FlorianD](https://discourse.julialang.org/u/FlorianD)\
**Post date:** [September 15, 2022, 5:45pm UTC](https://discourse.julialang.org/t/using-modelingtoolkit-jl-to-model-a-matrix-ode/82490/2 "2022-09-15T17:45:55Z")

</div>

I asked a similar question in this issue ["Axes of broadcast not known" with differential matrix equation · Issue #1740 · SciML/ModelingToolkit.jl · GitHub](https://github.com/SciML/ModelingToolkit.jl/issues/1740)

If your remove the `[]` around the equation and use `D.()` you can create the ODESystem.

```julia
using ModelingToolkit

using DifferentialEquations: solve

@variables x y[1:2, 1:2](x)

@parameters f g

mat = [1.0 2.0; 0.0 0.5]

D = Differential(x)

@named sys = ODESystem(D.(y) ~ (g * y + f * mat * y))

prob = ODEProblem(sys, [y => [1.0 0.0; 0.0 0.0]], [g => 1.0, f => 2.0])

```

But then your ODEProblem fails and I don’t know how to fix it.

```julia
ERROR: ArgumentError: SymbolicUtils.Symbolic[y[1, 1](x), y[2, 1](x), y[1, 2](x), y[2, 2](x)] are missing from the variable map.

```

---

<div class="post-metadata">

**Author:** ![FlorianD](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/floriand/32/38696_2.png) [@FlorianD](https://discourse.julialang.org/u/FlorianD)\
**Post date:** [September 15, 2022, 5:55pm UTC](https://discourse.julialang.org/t/using-modelingtoolkit-jl-to-model-a-matrix-ode/82490/3 "2022-09-15T17:55:59Z")

</div>

You also need to specify the timespan of your ODEProblem

This works, but its not exactly what you want. (I hard coded the values of parameters g and f in the equation)

```julia
using ModelingToolkit

using DifferentialEquations: solve

@variables x y[1:2, 1:2](x)

@parameters f g

mat = [1.0 2.0; 0.0 0.5]

D = Differential(x)

@named sys = ODESystem(D.(y) ~ (1.0 * y + 2.0 * mat * y))

t0 = 0.0

tf = 1.0

prob = ODEProblem(sys,[1.0 0.0; 0.0 0.0],(t0,tf))

solve(prob)

```

---

<div class="post-metadata">

**Author:** ![qwerty](https://avatars.discourse-cdn.com/v4/letter/q/4491bb/32.png) [@qwerty](https://discourse.julialang.org/u/qwerty)\
**Post date:** [September 17, 2022, 12:44pm UTC](https://discourse.julialang.org/t/using-modelingtoolkit-jl-to-model-a-matrix-ode/82490/4 "2022-09-17T12:44:04Z")

</div>

```julia
eqs = D.(y) ~ (g * y + f * mat * y)
@named sys = ODESystem(eqs)

y0 = collect(y .=> [1.0 0.0; 0.0 0.0])
p = [g => 1.0, f => 2.0]
tspan = (0.0, 0.5)

prob = ODEProblem(sys, y0, tspan, p)

```