  The exact factorization (EF) approach to coupled electron-ion dynamics
recasts the time-dependent molecular Schr\"odinger equation as two coupled
equations, one for the nuclear wavefunction and one for the conditional
electronic wavefunction. The potentials appearing in these equations have
provided insight into non-adiabatic processes, and new practical non-adiabatic
dynamics methods have been formulated starting from these equations. Here we
provide a first demonstration of a self-consistent solution of the exact
equations, with a preliminary analysis of their stability and convergence
properties. The equations have an unprecedented mathematical form, involving a
non-Hermitian Hamiltonian, and so the usual numerical methods for
time-dependent Schr\"odinger fail when applied in a straightforward way to the
EF equations. We find an approach that enables stable propagation long enough
to witness non-adiabatic behavior in a model system before non-trivial
instabilities take over. Implications for the development and analysis of
EF-based methods are discussed.
