In this paper, a rigorous method for the analysis of the propagation characteristics of various rectangular waveguides is proposed based on the recent development of the numerical mode matching method for multi-region, vertically stratified media. Unlike the conventional two-dimensional finite element method which is widely used for the analysis of dielectric waveguide problems, the numerical mode matching method proposed in this paper uses the propagation properties of the waves on the cross-section of the waveguides to convert the two-dimensional problem into a few one-dimensional problems. These one-dimensional problems are then solved with the one-dimensional finite element method, resulting in a great saving of the number of unknowns needed. A major difference between this method and the conventional mode matching method is that the eigenvalues and eigenmodes of the constituent layered media are found numerically by using a one-dimensional finite element method, which makes the method more versatile than the conventional mode matching method. This analysis method has been implemented numerically for analyzing arbitrary structures of rectangular dielectric waveguides. The comparison of the results with those available in literature are favorable. With the numerical results for various waveguide structures, the method is demonstrated to be efficient and capable of solving very complex structures such as coupled rib waveguides and coupled channel waveguides.