Guided and leaky modes of planar dielectric waveguides are eigensolutions of a singular Sturm-Liouville problem. This paper describes how this problem can be transformed into a quartic eigenvalue problem, which in turn can be converted into a generalized eigenvalue problem. Thus standard iterative methods, such as Arnoldi methods, can be used to compute the spectrum. We show how the shifts in the Arnoldi methods must be selected to obtain convergence to the dominant modes. In addition, by using high-order flnite elements, the resulting solutions can be made extremely accurate. Numerical examples demonstrate the speed and accuracy as well as the stability of the method. It is easy to characterize the modes of planar dielectric waveguides by an algebraic equation. There are numerous papers that discuss numerical methods that are based on flnding the roots of the characteristic function, either by using the argument principle of complex analysis (1{3,11), or by a continuation method (7). However, the numerical solution of the equation sufiers from numerical instabilities which are caused by the exponential scaling of the characteristic function. It is well known that for this reason the standard methods are often unreliable especially if there are many or thick layers present, if the frequency is large or if there are lossy layers, see, e.g., (5,9,12). In (12), we presented a new variational formulation of the Sturm-Liouville problem that char- acterizes the guided and leaky modes. After discretization, the variational form reduces to either a quadratic or a quartic matrix eigenvalue problem. Solving the eigenvalue problem is numerically stable, even for waveguides with thick layers or arbitrary index proflles. In (12), we used a low-order discretization scheme and solved the eigenvalue problem by a direct method. This demonstrated the general feasibility of the approach, but the convergence of the higher-order modes is slow and the high cost of the numerical linear algebra limits the applicability to relatively simple structures. In the present paper we address the slow convergence problem by using piecewise high-order polynomial flnite elements. Furthermore, we employ an iterative eigensolver that is capable to exploit the sparseness of the flnite-element matrices. To ensure that an iterative method converges to a speciflc eigenvalue one has to shift the problem such that the selected eigenvalue is dominant. We will present a strategy for selecting the shift such that all dominant guided and leaky modes are found. We will conclude with some numerical results that demonstrate the usefulness of our improved method. 2. VARIATIONAL FORMULATION In this section, we brie∞y review the derivation of the variational formulation that characterizes the modes of a dielectric waveguide. For more details we refer to (12). We consider an inflnite medium where the refractive index n(x) = p †(x) is a function of x and we let z be the direction of propagation. In this case the modes are of the form E(x;z;t) = exp(ii!t + iflz)`(x)^ (TE wave) or H(x;z;t) = exp(ii!t + iflz)`(x)^ (TM wave). To simplify notations we only discuss the TE waves, as the treatment of TM waves is completely analogous and can be found in (12). The lateral dependence of the mode is given by the functionwhich satisfles ` 00 (x) + i k 2 n 2 (x) i fl 2 ¢ `(x) = 0; x 2R;