In this paper we present solutions of evolution equations for inclusive distribution of gluons as produced by jet traversing quark–gluon plasma. We reformulate the original equations in such a form that virtual and unresolved-real emissions as well as unresolved collisions with medium are resummed in a Sudakov-type form factor. The resulting integral equations are then solved most efficiently with use of newly developed Markov Chain Monte Carlo algorithms implemented in a dedicated program called MINCAS. Their results for a gluon energy density are compared with an analytical solution and a differential numerical method. Some results for gluon transverse-momentum distributions are also presented. They exhibit interesting patterns not discussed so far in the literature, in particular a departure from the Gaussian behaviour – which does not happen in approximate analytical solutions.