  This paper shows how to use the shooting method, a classical numerical
algorithm for solving boundary value problems, to compute the Riemannian
distance on the Stiefel manifold $ \mathrm{St}(n,p) $, the set of $ n \times p
$ matrices with orthonormal columns. The proposed method is a shooting method
in the sense of the classical shooting methods for solving boundary value
problems; see, e.g., Stoer and Bulirsch, 1991. The main feature is that we
provide an approximate formula for the Fr\'{e}chet derivative of the geodesic
involved in our shooting method. Numerical experiments demonstrate the
algorithms' accuracy and performance. Comparisons with existing
state-of-the-art algorithms for solving the same problem show that our method
is competitive and even beats several algorithms in many cases.
