The modes of a planar multilayer waveguide can readily be determined by using an algorithm that exploits the analyticity of the waveguide dispersion relation. No prior knowledge of the number of solutions or their approximate locations is required; the algorithms finds all the solutions within the region of interest without risk of confusion because of degeneracy or near degeneracy. This method finds both leaky and bound modes in structures of transparent and lossy materials.