A rigorous modeling of vertical coaxial probes embedded in thick (uniaxial) anisotropic multilayered structures is presented. An integral equation for the current along the isolated probe is derived and solved by the method of moments. Two types of testing are used: point-matching and Galerkin. The steepest-descent method with a modified saddle-point approximation is used to obtain the asymptotic electromagnetic far field. Results are discussed in terms of the surface and leaky waves locations in the steepest-descent complex plane. The effect of anisotropy is pointed out. Computed results are compared with measurements.