Abstract This paper presents an accurate and efficient semianalytical method based on the Galerkin procedure for solving electromagnetic wave propagation problems in multilayer inhomogeneous cylindrical dielectric waveguides. The method represents the field in each inhomogeneous layer by a linear combination of eigenfunctions with unknown coefficients, which are expressed using the inner products of a series of basis functions, following the Galerkin procedure. The continuity of the field and its radial derivative is enforced at the interface between adjacent layers. By applying this procedure to all inhomogeneous layers, the Helmholtz equations are transformed into linear algebraic equations with expanded coefficients in matrix form, allowing the complicated wave propagation problem in a multilayer inhomogeneous waveguide to be solved as a matrix eigenvalue problem. The method is validated by providing detailed propagation characteristics for various multilayer inhomogeneous cylinders with different permittivity profiles. The accuracy and efficiency of the proposed method are demonstrated through comparisons with results obtained using other numerical techniques.