  A method for calculating phase-space master integrals for the decay process
$1 \to n$ massless partons in QCD using integration-by-parts and differential
equations techniques is discussed. The method is based on the appropriate
choice of the basis for master integrals which leads to significant
simplification of differential equations. We describe an algorithm how to
construct the desirable basis, so that the resulting system of differential
equations can be recursively solved in terms of (G)HPLs as a series in the
dimensional regulator $\epsilon$ to any order. We demonstrate its power by
calculating master integrals for the NLO time-like splitting functions and
discuss future applications of the proposed method at the NNLO precision.
