Summary: Many applications give rise to dynamical systems in the form of a vector field with a phase space of moderate dimension. Examples are the Lorenz equations, mechanical and other oscillators, and models of spiking neurons. The global dynamics of such a system is organized by the stable and unstable manifolds of the saddle points, of the saddle periodic orbits, and, more generally, of all compact invariant manifolds of saddle type. Except in very special circumstances, the (un)stable manifolds are global objects that cannot be found analytically but need to be computed numerically. This is a nontrivial task when the dimension of the manifold is larger than one. Here, we present an algorithm to compute the \(k\)-dimensional unstable manifold of an equilibrium or periodic orbit (or a more general normally hyperbolic invariant manifold) of a vector field with an \(n\)-dimensional phase space\, with \(1< k < n\). Stable manifolds are computed by considering the flow for negative time. The key idea is to view the unstable manifold as a purely geometric object, hence disregarding the dynamics on the manifold, and compute it as a list of approximate geodesic level sets, which are (topological) (\(k-1\))-spheres. Starting from a (\(k-1\))-sphere in the linear eigenspace of the equilibrium or periodic orbit, the next geodesic level set is found in a local (and changing) coordinate system given by hyperplanes perpendicular to the last geodesic level set. In this setup, the mesh points defining the approximation of the next geodesic level set can be found by solving boundary value problems. By appropriately adding or removing mesh points it is ensured that the mesh that represents the computed manifold is of a prescribed quality. The algorithm is presently implemented to compute two-dimensional manifolds in a phase space of arbitrary dimension. In this case, the geodesic level sets are topological circles and the manifold is represented as a list of bands between consecutive level sets. We use color to distinguish between consecutive bands or to indicate geodesic distance from the equilibrium or periodic orbit, and we also show how geodesic level sets change with increasing geodesic distance. This is very helpful when one wants to understand the often very complicated embeddings of two-dimensional (un)stable manifolds in phase space. The properties and performance of our method are illustrated with several examples, including the stable manifold of the origin of the Lorenz system, a two-dimensional stable manifold in a four-dimensional phase space arising in a problem in optimal control, and a stable manifold of a periodic orbit that is a Möbius strip. Each illustration is accompanied by an animation (supplied with this paper).