We propose to compute approximations to invariant sets in dynamical systems by minimizing an appropriate distance between a suitably selected finite set of points and its image under the dynamics. We demonstrate, through computational experiments, that this approach can successfully converge to approximations of (maximal) invariant sets of arbitrary topology, dimension, and stability, such as, e.g., saddle type invariant sets with complicated dynamics. We further propose to extend this approach by adding a Lennard-Jones type potential term to the objective function, which yields more evenly distributed approximating finite point sets, and illustrate the procedure through corresponding numerical experiments. In the phase space of any nonlinear dynamical system, the “skeleton” of the global dynamical behavior consists of the invariant sets of the system, e.g., fixed points, periodic orbits, general recurrent sets, and the connecting orbits/invariant manifolds between them. Computational methods for approximating invariant sets have been, and will continue to be, a major part of the “toolkit” of every dynamical systems researcher, whether on the mathematical or on the modeling side. In this contribution we devise and implement a new variational approach for this task, which is able to compute invariant sets of arbitrary dimension, topology, and stability type. In addition—and in contrast to classical techniques—our method provides an approximate parametrization of the invariant set, which can be (smoothly) followed in parameter space. I. INTRODUCTION