Von Neumann Equation Density Matrix Implementation

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