  In this paper, we propose a numerical method to solve the classic
$L^2$-optimal transport problem. Our algorithm is based on use of multiple
shooting, in combination with a continuation procedure, to solve the boundary
value problem associated to the transport problem. We exploit the viewpoint of
Wasserstein Hamiltonian flow with initial and target densities, and our method
is designed to retain the underlying Hamiltonian structure. Several numerical
examples are presented to illustrate the performance of the method.
