We propose a new numerical method for exactly determining modes of a general multilayer waveguide. The proposed method has two features: simple implementation and fast calculation. Its implementation is simple since it is based on a complex root-finding algorithm reported recently. This algorithm is very intuitive and consists of the process of isolating each zero and that of determining the isolated zeros. To reduce calculation time, we propose and use an efficient scheme for evaluations of the eigenvalue function and its derivative. Using this scheme, we can decrease the number of arithmetic operations necessary for the analysis compared to the existing scheme based on the transfer matrix method. The proposed method is applied to analyze lossless and lossy waveguides and an anti-resonant reflective optical waveguide. The analysis results are compared to those in literatures and the comparisons show good agreement. In addition, it is shown that the proposed method is useful in analyzing the effects of a metal layer or a high index substrate on the integrated optical devices.