The bivariate high-temperature expansion of the spin-spin correlation function of the three-dimensional classical XY ͑planar rotator͒ model, with spatially anisotropic nearest-neighbor couplings, is extended from the tenth through the 21st order. The computation is carried out for the simple cubic lattice, in the absence of magnetic field, in the case in which the coupling strength along the z axis of the lattice is different from those along the x and the y axes. It is then possible to determine accurately the critical temperature as function of the parameter R which characterizes the coupling anisotropy and to check numerically the universality, with respect to R, of the critical exponents of the three-dimensional anisotropic system. The analysis of our data also shows that the main predictions of the generalized scaling theory for the crossover from the three-dimensional to the two-dimensional critical behavior are compatible with the series extrapolations.