We extend the auxiliary-mass-flow (AMF) method originally developed for Feynman loop integration to calculate integrals involving also phase-space integration. Flow of the auxiliary mass from the boundary ($infty$) to the physical point ($0^+$) is obtained by numerically solving differential equations with respective to the auxiliary mass. For problems with two or more kinematical invariants, the AMF method can be combined with traditional differential equation method by providing systematical boundary conditions and highly nontrivial self-consistent check. The method is described in detail with a pedagogical example of $e^+e^-rightarrow gamma^* rightarrow tbar{t}+X$ at NNLO. We show that the AMF method can systematically and efficiently calculate integrals to high precision.