Markov-modulated Brownian motion is a popular tool to model continuous-time phenomena in a stochastic context. The main quantity of interest is the invariant density, which satisfies a differential equation associated with the quadratic matrix polynomial $P(z) = Vz^2-Dz +Q$, where the matrices $V$ and $D$ are diagonal and $Q$ is the transition matrix of a discrete-time Markov chain. Its solution is typically constructed by computing an invariant pair of $P(z)$ associated with its eigenvalues in the left half-plane, or by solving the matrix equation $X^2V-XD+Q=0$. We show that these tasks can be solved using a componentwise accurate algorithm based on Cyclic Reduction, generalizing the recently appeared algorithms for the linear case ($V=0$). We give a proof of the numerical stability of our algorithm in the componentwise sense; the same proof applies to Cyclic Reduction in a more general M-matrix setting which appears in other applications such as the modelling of QBD processes.