Internal flows inside gravitationally stable astrophysical objects, such as the Sun, stars and compact stars are compressed and extremely subsonic. Such low Mach number flows are usually encountered when studying for example dynamo action in stars, planets, the hydro-thermodynamics of X-ray bursts on neutron stars and dwarf novae. Treating such flows is numerically complicated and challenging task. We aim to present a robust numerical tool that enables modeling the time-evolution or quasi-stationary of stratified low Mach number flows under astrophysical conditions. It is argued that astrophysical low Mach number flows cannot be considered as an asymptotic limit of incompressible flows, but rather as highly compressed flows with extremely stiff pressure terms. Unlike the pseudo-pressure in incompressible fluids, a Possion-like treatment for the pressure would smooth unnecessarily the physically induced acoustic perturbations, thereby violating the conservation character of the compressible equations. Moreover, classical dimensional splitting techniques, such as ADI or Line-Gauss-Seidel methods are found to be unsuited for modeling compressible flows with low Mach numbers. In this paper we present a nonlinear Newton-type solver that is based on the defect-correction iteration procedure and in which the Approximate Factorization Method (AFM) is used as a preconditioner. This solver is found to be sufficiently robust and is capable of capturing stationary solutions for viscous rotating flows with Mach number as small as $\mcal{M} \approx 10^{-3},$ i.e., near the incompressibility limit.