This set of coupled equations may thus be propagated using an ODE integrator. The time integration is performed using the well-known adaptive time-step Bulirsch-Stoer scheme with error tolerance set to 10 -9 . Moreover, an upper bound is set on the integration step so that no neighboring trajectories spacing should decrease by a factor above 40% from one time step to the next. For comparison, we also performed a wavefunction-based calculation of the time evolution of the corresponding TDSE, using the Crank-Nicholson algorithm.