We review the numerical techniques for ideal and non-ideal magneto-hydrodynamics (MHD) used in the context of star formation simulations. We outline the specific challenges offered by modelling star forming environments, which are dominated by supersonic turbulence in a radiative, self-gravitating fluid. These conditions are rather unique in physics and engineering and pose particularly severe restrictions on the robustness and accuracy of numerical codes. One striking aspect is the formation of collapsing fluid elements leading to the formation of singularities that represent point-like objects, namely the proto-stars. Although a few studies have attempted to resolve the formation of the first and second Larson's cores, resolution limitations force us to use sink particle techniques, with sub-grid models to compute the accretion rates of mass, momentum and energy, as well as their ejection rate due to radiation and jets from the proto-stars. We discuss the most popular discretisation techniques used in the community, namely smoothed particle hydrodynamics, finite difference and finite volume methods, stressing the importance to maintain a divergence free magnetic field. We discuss how to estimate the truncation error of a given numerical scheme, and its importance in setting the magnitude of the numerical diffusion. This can have a strong impact on the outcome of these simulations, where viscosity and resistivity are implemented at the grid scale. We then present various numerical techniques to model non-ideal MHD effects (Ohmic and ambipolar diffusion, and Hall effect). These physical ingredients are posing strong challenges in term of resolution and time stepping. For the latter, several strategies are discussed to overcome the limitations due to ridiculously small time steps. An important aspect of star formation simulations is the radiation field. We discuss the current state-of-the-art, with a variety of techniques offering pros and cons in different conditions. Finally, we present more advanced strategies to mitigate the adverse effect of finite numerical resolution, which are very popular in the context of supersonic, self-gravitating fluids, namely adaptive mesh refinement, moving meshes, SPH and high-order methods. Advances in these directions are likely to trigger immense progresses in the future of our field.