Many devices and structures used to guide electromagnetic waves are conformal with the cylindrical coordinates. Sensitive applications of microwave engineering and integrated optical devices often use non-homogeneous, anisotropic and dissipative materials, so that the research for robust and accurate computational models is a topic of remarkable interest for Electrical Engineering. This work presents a semi-analytical technique for solving boundary-value problems associated with cylindrical, anisotropic, and non-homogeneous waveguides. Our methodology allows us to model structures with radial layers, with uniaxial anisotropy, and with losses. The proposed solution starts from Maxwell s equations for time-harmonic electromagnetic fields and employs a modal expansion in terms of the Bessel-Fourier series. The eigenvalues associated with the problem are obtained using the winding number method, in which several approaches for calculating complex-plane contour integrals are explored in detail. In order to properly analyze the junctions between sections of stratified waveguides, we employ a mode-matching technique based on the conservation of the Reaction of the fields. Our formulation can handle the effects of excitation and coupling between pure modes (TM, TE, and TEM) in homogeneous waveguides, as well as hybrid modes in complex structures. A series of numerical results are presented and show the capacity of the methodology developed here to correctly characterize cylindrical structures composed of complex media (inhomogeneous, anisotropic, and dissipative) in a robust and computationally-efficient fashion if compared to other conventional computational electromagnetic techniques.