In this work, we introduce a couple of algorithms to compute the stationary probability distribution for the chemical master equation (CME) of arbitrary chemical networks. We further find the conditions guaranteeing the algorithms' convergence and the unity and stability of the stationary distribution. Next, we employ these algorithms to study the mRNA and protein probability distributions in a gene regulatory network subject to negative feedback regulation. In particular, we analyze the influence of the promoter activation/deactivation speed on the shape of such distributions. We find that a reduction of the promoter activation/deactivation speed modifies the shape of those distributions in a way consistent with the phenomenon known as mRNA (or transcription) bursting.