I’m trying to implement the von Neumann Equation for a given 4×4 density Matrix with a time dependent Hamiltonian `Hp[t_]`

in Mathematica but I get stuck.

`Format[y[a__]] := Subscript[y, a] rho[t_] := Array[x[##][t] &, {4, 4}] sol = NDSolve[{I*rho'[t] == Hp[t].rho[t] - rho[t].Hp[t], rho[0] == rhoIni}, {rho}, {t, 0, 10}] `

However I only get the output

`{{rho -> rho}} `

So I guess something is structurally wrong with my code. I try to extract a solution by writing

`rho[t_] = rho[t] /. sol `

But this doesn’t work as there is no solution anyways. Maybe you can help me

Thanks in advance