While magnetic fields have long been considered significant for the evolution of magnetic non-degenerate stars and compact stars, it has become clear in recent years that, in fact, all stars are deeply affected by their effects. This is particularly true regarding their internal angular momentum distribution, but magnetic fields may also influence internal mixing processes and even the fate of the star. We propose a new framework for stellar evolution simulations in which the interplay between magnetic field, rotation, mass loss, and changes in the stellar density and temperature distributions are treated self-consistently. For average large-scale stellar magnetic fields that are symmetric to the axis of the rotation of the star, we derive 1D evolution equations for the toroidal and poloidal components from the mean-field magnetohydrodynamic equation by applying Alfvén’s theorem; and, hence, a conservative form of the angular momentum transfer due to the Lorentz force is formulated. We implement our formalism into a numerical stellar evolution code and simulate the magneto-rotational evolution of 1.5 M⊙ stars. The Lorentz force aided by the Ω effect imposes torsional Alfvén waves propagating through the magnetized medium, leading to near-rigid rotation within the Alfvén timescale. Our models, with different initial spins and B-fields, can reproduce the main observed properties of Ap/Bp stars. Calculations that are extended to the red-giant regime show a pronounced core-envelope coupling, which are capable of reproducing the core and surface rotation periods already determined by asteroseismic observations.