Crystal and continuum plasticity theories describe plastic behavior in metals on two tightly linked material scales, i.e. the grain scale and the component scale. From the multi-scale material modelling perspective, the efficiency and robustness of the solvers involved in the modelling framework, are of vital importance for the overall performance. Stable, fully implicit return-mapping algorithms are developed for solving both crystal and continuum plasticity models. Their implementation into user-defined material subroutines (UMAT) in a finite-element (FE) software are made available as open-source. Numerical stability is gained by an improved initial guess for the stress solution and by applying a line search for each Newton iteration. FE simulations were run, demonstrating the performance of the new implementations. A simulation of the necking of an aluminium single crystal and of the deformation of an aluminium polycrystalline representative volume element were performed, demonstrating the stability and high efficiency of the new crystal plasticity approach. Similarly, continuum plasticity simulations revealed unconditional stability and very high efficiency of the new implementation of the advanced anisotropic Yld2004-18p yield function. It performs equally fast as the much simpler von Mises and Hill implementations available as standard models in the software. This enables the full exploitation of advanced yield functions in industrial FE applications and sets the new standard for the metal forming industry.