Summary: Computation of electromagnetic fields due to point sources (Hertzian dipoles) in cylindrically stratified media is a classical problem for which analytical expressions of the associated tensor Green’s function have been long known. However, under finite-precision arithmetic, direct numerical computations based on the application of such analytical (canonical) expressions invariably lead to underflow and overflow problems related to the poor scaling of the eigenfunctions (cylindrical Bessel and Hankel functions) for extreme arguments and/or high-order, as well as convergence problems related to the numerical integration over the spectral wavenumber and to the truncation of the infinite series over the azimuth mode number. These problems are exacerbated when a disparate range of values is to be considered for the layers’ thicknesses and material properties (resistivities, permittivities, and permeabilities), the transverse and longitudinal distances between source and observation points, as well as the source frequency. To overcome these challenges in a systematic fashion, we introduce herein different sets of range-conditioned, modified cylindrical functions (in lieu of standard cylindrical eigenfunctions), each associated with nonoverlapped subdomains of (numerical) evaluation to allow for stable computations under any range of physical parameters. In addition, adaptively-chosen integration contours are employed in the complex spectral wavenumber plane to ensure convergent numerical integration in all cases. We illustrate the application of the algorithm to problems of geophysical interest involving layer resistivities ranging from \(1000 \Omega \text{m}\) to \(10^{- 8} \Omega \text{m}\), frequencies of operation ranging from 10 MHz down to the low magnetotelluric range of 0.01 Hz, and for various combinations of layer thicknesses.