A few years ago I had a similar problem regarding a system of complex equations and I had found a workaround: ModelingToolkit matrix form system of equations and complex values.

As far as I can tell complex symbolics are still not supported in ModelingToolkit, but now unfortunately my workaround is no longer functioning either. Some of the ModelingToolkit syntax has changed in the almost 3 years since, so here is an updated MWE:

```
@parameters icomplex Γ32 δl δμ Ωl Ωμ
@variables t (ρ(t))[1:3,1:3]
D = Differential(t)
H = [
-δl-δμ -Ωμ/2 0
-Ωμ/2 -δl -Ωl/2
0 -Ωl/2 0
]
L = [
0 0 -Γ32*ρ[1,3]/2
0 Γ32*ρ[3,3] -Γ32*ρ[2,3]/2
-Γ32*ρ[3,1]/2 -Γ32*ρ[3,2]/2 -Γ32*ρ[3,3]
]
eqns = D(ρ) ~ -icomplex*(H*ρ-ρ*H) + L
@mtkbuild bloch = ODESystem(Symbolics.scalarize(eqns), t)
ρᵢ = zeros(ComplexF64,3,3)
ρᵢ[2,2] = 1
u0 = [ρ[idx,idy] => ρᵢ[idx,idy] for idx in 1:3 for idy in 1:3]
p = [Ωl => 1, Ωμ=> 0, δl => 0., δμ => 0., Γ32 => 1., icomplex => 1.0im]
prob = ODEProblem(bloch, u0, (0., 20.), p, jac = true)
sol = solve(prob, Tsit5())
```

The issue occurs when constructing the ODEProblem, throwing the an `InexactError: Float64(0.0 + 1.0im)`

. It seems like internally it tries to convert `icomplex`

to a `Float64`

.

Is there a different workaround possible?