One of the most accurate methods for solving the time-dependent Schrödinger equation uses a combination of the dynamic Fourier method with the split-operator algorithm on a tensor-product grid. To reduce the number of required grid points, we let the grid move together with the wavepacket but find that the naïve algorithm based on an alternate evolution of the wavefunction and grid destroys the time reversibility of the exact evolution. Yet, we show that the time reversibility is recovered if the wavefunction and grid are evolved simultaneously during each kinetic or potential step; this is achieved by using the Ehrenfest theorem together with the splitting method. The proposed algorithm is conditionally stable, symmetric, and time-reversible and conserves the norm of the wavefunction. The preservation of these geometric properties is shown analytically and demonstrated numerically on a three-dimensional harmonic model and collinear model of He–H2 scattering. We also show that the proposed algorithm can be symmetrically composed to obtain time-reversible integrators of an arbitrary even order. We observed 10 000-fold speedup by using the tenth-order instead of the second-order method to obtain a solution with a time discretization error below 10−9. Moreover, using the adaptive grid instead of the fixed grid resulted in a 64-fold reduction in the required number of grid points in the harmonic system and made it possible to simulate the He–H2 scattering for six times longer while maintaining reasonable accuracy. Applicability of the algorithm to high-dimensional quantum dynamics is demonstrated using the strongly anharmonic eight-dimensional Hénon–Heiles model.