A computational scheme, based on a time-dependent extension of non-collinear spin density functional theory, for the simultaneous simulation of charge and magnetization dynamics in molecular systems is presented. We employ a second-order Magnus propagator combined with an efficient predictor-corrector scheme that allows us to treat large molecular systems over long simulation periods. The method is benchmarked against the low-frequecy dynamics of the H-He-H molecule where the magnetization dynamics can be modeled by the simple classical magnetization precession of a Heisenberg-Dirac-van Vleck Hamiltonian. Furthermore, the magnetic exchange couplings of the bimetallic complex [Cu(bpy)(H₂O)(NO₃)₂(μ-C₂O₄)] (BISDOW) are extracted from the low-frequency precession dynamics showing good agreement with coupling obtained from ground state energy differences. Our approach opens the possibility to perform real-time simulation of spin-related phenomena using time-dependent density functional theory in realistic molecular systems.