The guided and leaky modes of a planar dielectric waveguide are eigensolutions of a singular Sturm-Liouville problem. The modes are the roots of a characteristic function which can be found with several methods that have been introduced in the past. However, the evaluation of the characteristic function suffers from numerical instabilities, and hence it is often difficult to find all modes in a given range. Here a new variational formulation is introduced, which, after discretization, leads either to a quadratic or a quartic eigenvalue problem. The modes can be computed with standard software for polynomial eigenproblems. Numerical examples show that the method is numerically stable and guarantees a complete set of solutions.