  Dynamics of an open $N$-state quantum system is typically modeled with a
Markovian master equation describing the evolution of the system's density
operator. By using generators of $SU(N)$ group as a basis, the density operator
can be transformed into a real-valued 'Bloch vector'. The Lindbladian, a
super-operator which serves a generator of the evolution, %in the master
equation, can be expanded over the same basis and recast in the form of a real
matrix. Together, these expansions result is a non-homogeneous system of
$N^2-1$ real-valued linear differential equations for the Bloch vector. Now one
can, e.g., implement a high-performance parallel simplex algorithm to find a
solution of this system which guarantees exact preservation of the norm and
Hermiticity of the density matrix. However, when performed in a straightforward
way, the expansion turns to be an operation of the time complexity
$\mathcal{O}(N^{10})$. The complexity can be reduced when the number of
dissipative operators is independent of $N$, which is often the case for
physically meaningful models. Here we present an algorithm to transform quantum
master equation into a system of real-valued differential equations and
propagate it forward in time. By using a scalable model, we evaluate
computational efficiency of the algorithm and demonstrate that it is possible
to handle the model system with $N = 10^3$ states on a single node of a
computer cluster.
