  The explicit split-operator algorithm is often used for solving the linear
and nonlinear time-dependent Schr\"{o}dinger equations. However, when applied
to certain nonlinear time-dependent Schr\"{o}dinger equations, this algorithm
loses time reversibility and second-order accuracy, which makes it very
inefficient. Here, we propose to overcome the limitations of the explicit
split-operator algorithm by abandoning its explicit nature. We describe a
family of high-order implicit split-operator algorithms that are
norm-conserving, time-reversible, and very efficient. The geometric properties
of the integrators are proven analytically and demonstrated numerically on the
local control of a two-dimensional model of retinal. Although they are only
applicable to separable Hamiltonians, the implicit split-operator algorithms
are, in this setting, more efficient than the recently proposed integrators
based on the implicit midpoint method.
