  This work is a continuation of our efforts to develop an efficient implicit
solver for multidimensional hydrodynamics for the purpose of studying important
physical processes in stellar interiors, such as turbulent convection and
overshooting. We present an implicit solver that results from the combination
of a Jacobian-Free Newton-Krylov method and a preconditioning technique
tailored to the inviscid, compressible equations of stellar hydrodynamics. We
assess the accuracy and performance of the solver for both 2D and 3D problems
for Mach numbers down to $10^{-6}$. Although our applications concern flows in
stellar interiors, the method can be applied to general advection and/or
diffusion-dominated flows. The method presented in this paper opens up new
avenues in 3D modeling of realistic stellar interiors allowing the study of
important problems in stellar structure and evolution.
