  Impurities coupled to superconductors offer a controlled platform to
understand the interplay between superconductivity, many-body interactions, and
non-equilibrium physics. In the equilibrium situation, local interactions at
the impurity induce a transition between the spin-singlet to the spin-doublet
ground state, resulting in a supercurrent sign reversal ($0-\pi$ transition).
In this work, we apply the exact time-dependent density matrix renormalization
group method to simulate the transient dynamics of such superconducting
systems. We also use a perturbative approximation to analyze their properties
at longer times. These two methods agree for a wide range of parameters. In a
phase-biased situation, the system gets trapped in a metastable state
characterized by a lower supercurrent compared to the equilibrium case. We show
that local Coulomb interactions do not provide an effective relaxation
mechanism for the initially trapped quasiparticles. In contrast, other
relaxation mechanisms, such as coupling to a third normal lead, make the
impurity spin relax for parameter values corresponding to the equilibrium $0$
phase. For parameters corresponding to the equilibrium $\pi$ phase the impurity
converges to a spin-polarized stationary state. Similar qualitative behavior is
found for a voltage-biased junction, which provides an effective relaxation
mechanism for the trapped quasiparticles in the junction.
