I am trying to use the ModelingToolkit together with Unitful to model a coupled system of two ODEs with different units. The following code fails on ODEProblem with a DimensionError unless both units in u0 are the same:

```
using ModelingToolkit
using DifferentialEquations
using Plots
using Unitful
@parameters t k1 k2
@variables y(t) z(t)
@derivatives D'~t
eqs = [D(y) ~ k1*y, D(z) ~ k2*z]
de = ODESystem(eqs)
u0 = [y => 1.5u"V", z => 3.1u"N"]
p = [k1 => 1.0u"1/s", k2 => 1.0u"1/s"]
tspan = (0.0u"s",1e-5u"s")
prob = ODEProblem(de,u0,tspan, p, jac=false) # ERROR: DimensionError: V and 3.1 N are not dimensionally compatible.
sol = solve(prob, Euler(), dt=1e-6u"s")
```

I know that heterogeneous variables are an issue, as discussed here: DifferentialEquations.jl and systems with heterogenous units

The problem in that case was the adaptive stepping that can’t calculate a vector norm if the vector has heterogeneous units.

The suggested solution was to use ArrayPartition:

ODE systems with heterogenous units · Issue #147 · SciML/DifferentialEquations.jl · GitHub

It suspect that I am running into a similar problem, even though in my case the code fails at ODEProblem(), not on solve().

What is the best approach to use units in my case?