Summary: A new numerical approach based on consistent operator splitting is presented for computing compressible, highly stratified flows in astrophysics. The algorithm is particularly designed to search for steady or almost steady solutions for the time-dependent Navier-Stokes equations, describing viscous flow under the influence of a strong gravitational field. The algorithm proposed is multidimensional and works in Cartesian, cylindrical or spherical co-ordinates. It uses a second-order finite volume scheme with third-order upwinding and a second-order time discretization. We use an adaptive time step control and monotonic multilevel grid distribution to speed up convergence. This method has been incorporated into a hydrodynamical code by which, for the first time, for two-dimensional models the dynamics of the boundary layer in the accretion disk around a compact star could be computed over the whole viscous time scale.