We have formulated an efficient algorithm for propagating excited-state wavefunctions involving many electronic states using well-established overlap-based methods and phaseconventions. In particular, we have combined an exact method for constructing the overlap matrix for CIS wavefunctions using a bi-orthogonal basis with an optimization protocol for choosing the proper phase of adiabatic states. We applied the resulting algorithm to investigate the high-lying excited-state dynamics of a large electronic system with many trivial crossings and degenerate states. Our results show that use of the bi-orthogonal basis can The final spectrum at time t = 145 fs, significantly after the scattering event. Note that the strongest peak (≈ 3.41 eV) is no longer a collective excitation; the plasmon has moved to 3.57 eV. accelerate the construction of the overlap matrix by orders of magnitude at each time step, and that choosing the optimal phase of the overlap matrix allows can allow for a larger classical time step for the propagation. Altogether, large savings can be accrued.