  One of the most accurate methods for solving the time-dependent
Schr\"{o}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\"ive 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, 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-H$_{2}$ scattering.
We also show that the proposed algorithm can be symmetrically composed to
obtain time-reversible integrators of an arbitrary even order. We observed
$10000$-fold speedup by using the tenth- 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-H$_{2}$ 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\'{e}non--Heiles model.
