Electromagnetic wave propagation problems related to multilayer inhomogeneous cylindrical dielectric waveguides are solved by using an accurate and efficient semianalytical method based on the Galerkin procedure. Each inhomogeneous layer's field is represented as a linear combination of eigenfunctions with unknown coefficients, and the eigenfunctions are expressed by the inner products of a series of basis functions in term of the Galerkin procedure. The continuity of the field and its radial derivative is applied at the interface between the adjacent layers. Repeating the same procedure for all inhomogeneous layers, the Helmholtz equations are transformed into linear algebraic equations of expanded coefficients in matrix form. Therefore, the complicated problem of wave propagation in a multilayer inhomogeneous waveguide could be reduced to a simple problem of solving a matrix's eigenvalues. Detailed propagation characteristics are given for several representative and practical multilayer inhomogeneous cylinders with different kinds of permittivity profiles. The accuracy and the efficiency of the proposed method is then validated by comparing the calculated results with those obtained by using other numerical techniques.