I’ve started looking into ModelingToolkit and I wanted to solve optical bloch differential equations with it, which have a complex part.

I’ll preface this by saying I don’t have much experience with Julia yet, I started working with it recently.

A simple example of such a system would be a three level system with a hamiltonian

H = \begin{pmatrix}
-\delta_l - \delta_\mu & -\Omega_\mu/2 & 0 \\
\bar\Omega_\mu/2 & -\delta_l & -\Omega_l/2 \\
0 & -\bar\Omega_l/2 & 0
\end{pmatrix}

and a dissipator

L = \begin{pmatrix}
0 & 0 & -\Gamma_{32} \rho_{1,3}(t)/2 \\
0 & \Gamma_{32}\rho_{3,3}(t) & -\Gamma_{32}\rho_{2,3}(t)/2 \\
-\Gamma_{32}\rho_{3,1}(t)/2 & -\Gamma_{32}\rho_{3,2} & -\Gamma_{32}\rho_{3,3}(t)
\end{pmatrix}

The governing equation is:

\frac{\partial \rho(t)}{\partial t} = -i[H,\rho] + L

I can define all the variables and parameters in Julia and make the matrices

```
using ModelingToolkit
@variables t ρ[1:3,1:3](t)
@parameters Ωl Ωμ δl δμ Γ32
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]]
```

I run into trouble when setting up the `ODESystem`

in matrix notation.

```
eq = -1im*(H*ρ - ρ*H) + L
bloch = ODESystem(D(ρ) ~ eq)
```

This throws a conversion error (`Cannot convert an object of type Missing to Sym`

).

I then tried just writing it as a list of equations:

```
eqns = [D(ρ[idx,idy]) ~ eq[idx,idy] for idx in 1:3 for idy in 1:3]
bloch = ODESystem(eqns)
```

This however doesn’t work because of the complex part of the system, `eqns`

is an array of arrays, where the 1st entry is the real part and the 2nd entry the complex part, which then obviously doesn’t work with `ODESystem`

.

Finally as a test, when I use as an equation \frac{\partial \rho(t)}{\partial t} = -[H,\rho] + L (so no imaginary part) and use the list of equations approach I can create an `ODESystem`

and solve the ODE with `DifferentialEquations.jl`

, but it obviously doesn’t return any physical results then.

As a workaround I can use `build_function`

on `H`

and `L`

and manually create the rhs function for the solver, but I’d like to use `ModelingToolkit.jl`

because more complex examples will have \delta_l or other parameters be time dependent (so switched to variables), and this is rather painful to do when using the build_function approach because it requires manually creating the rhs function for each change.

To summarize, I have two main questions:

- Can you give
`ODESystem`

a matrix form of a system of equations - Is it possible to work with complex values for an
`ODESystem`

EDIT:

For completeness I’ll post the workaround I found (without matrix form or complex symbols). It basically uses another parameter to specify the complex part, which works with this simple example when the rates aren’t complex.

```
@parameter icomplex
eq = -icomplex*Symbolics.scalarize(H*ρ-ρ*H) + L
eqns = [D(ρ[idx,idy]) ~ eq[idx,idy] for idx in 1:3 for idy in 1:3]
bloch = ODESystem(eqns)
ρᵢ = 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 => 1im]
```

This can then be used with `DifferentialEquations.jl`

```
prob = ODEProblem(bloch, u0, (0., 200.), p, jac = true)
sol = solve(prob, Tsit5())
```