Multidimensional simulations of magnetohydrodynamic phenomena occurring in stellar interiors are essential for understanding how stars evolve and die. The highly subsonic flow regimes found in the regions deep inside stars pose severe challenges to conventional methods of computational MHD, such as the popular "high-resolution shock-capturing'' schemes. After giving a brief overview of work on astrophysical simulations (including also supernova explosions and common-envelope evolution) in our group at Heidelberg, we summarize the challenges and present suitable numerical solvers optimized for magnetized, low-Mach-number stellar flows, implemented in our Seven-League Hydro code. We show how the choice of the numerical method can drastically affect both the performance of the code and its accuracy in real astrophysical simulations.