  We propose a computationally efficient method to solve the dynamics of
operators of bosonic quantum systems coupled to their environments. The method
maps the operator under interest to a set of complex-valued functions, and its
adjoint master equation to a set of partial differential equations for these
functions, which are subsequently solved numerically. In the limit of weak
coupling to the environment, the mapping of the operator enables storing the
operator efficiently during the simulation, leading to approximately quadratic
improvement in the memory consumption compared with the direct approach of
solving the adjoint master equation in the number basis, while retaining the
computation time comparable. Moreover, the method enables efficient
parallelization which allows to optimize for the actual computational time to
reach an approximately quadratic speed up, while retaining the memory
consumption comparable to the direct approach. We foresee the method to prove
useful, e.g., for the verification of the operation of superconducting quantum
processors.
