A spectral-element type approach to determine the propagation characteristics of the bound modes of an open planar layered anisotropic waveguide is proposed. The main feature is the expansion of the electromagnetic field in each internal layer on different sets of Legendre polynomials. In the two embedding halfspaces, sets of weighted Laguerre polynomials are used. In this way, the method converges with exponential rate with increasing number of basis functions. No transcendental equation has to be solved and even modes that may exhibit coinciding or very close propagation constants are computed with great accuracy. A detailed analysis of the convergence properties of the algorithm is carried out and a way to estimate and control the error on the propagation constants is discussed.