An efficient and systematic strategy to find the propagation constants of real, complex, and leaky modes of covered and uncovered planar multilayered isotropic/uniaxial waveguides is presented. This strategy first builds up a pole-free characteristic function and then makes use of a root-searching scheme based on an integral nature method to search for its zeros. For uncovered waveguides, the branch points of the characteristic function can be removed by introducing the upper half-space vertical wavenumber as the working variable. Examples of the efficiency, reliability, and robustness of the presented technique are given for both covered and uncovered waveguides, which sets up this technique as a very convenient CAD tool. The method is also applied to study a new and interesting problem: the evolution of the modes of a grounded dielectric waveguide when increasing the permittivity of the upper semi-infinite medium.