We extend the auxiliary-mass-flow(AMF)method originally developed for Feynman loop integration to calculate integrals which also involve phase-space integration.The flow of the auxiliary mass from the boundary(oo)to the physical point(0+)is obtained by numerically solving differential equations with respective to the auxili-ary mass.For problems with two or more kinematical invariants,the AMF method can be combined with the tradi-tional differential-equation method,providing systematic boundary conditions and a highly nontrivial self-consist-ency check.The method is described in detail using a pedagogical example of e+e-→ y* →tt+X at NNLO.We show that the AMF method can systematically and efficiently calculate integrals to high precision.