A numerical scheme is suggested for accurate large-scale quantum dynamics simulations. The time-dependent Schrodinger equation with finite time-dependent interaction terms is replaced by an inhomogeneous equation with imaginary boundary operators applied along the time axis. This equation is solved globally for a finite time interval using recent Krylov subspace-based iterative methods that are accelerated by a Fourier grid preconditioner. The same scheme is applied also to time-independent reactive-scattering calculations with absorbing boundary operators where the operation of the Green’s function is carried out by solving an inhomogeneous time-independent equation. The scheme is economic in terms of both memory requirement and computation time, It is especially favorable when high grid densities are required, e.g., for representation of highly oscillatory fields or high-energy wave functions. Illustrative applications are given for representative models of bound and dissociative systems driven by time-dependent pulsed fields, and for time-independent calculations of the cumulative reaction probability for the generic reaction H+H-2 at high collision energies.