We have performed extensive Monte Carlo (MC) simulations in the NVT ensemble of the pair correlation function (p.c.f.) g(r) for square-well (SW) fluids for several densities, temperatures and potential widths. In addition, using the procedure developed by Smith et al. [J. Chem. Phys. 55 (1971) 4027], we have determined by MC-NVT simulations the zero- and first-order terms in the expansion of the p.c.f. in power series of the inverse of the reduced temperature at the same densities. Furthermore, we have compared the values of g(r) obtained from the expansion truncated beyond the first-order term, with those obtained directly by MC simulations performed on the SW system. The aim was mainly to test whether the perturbation expansion converges quickly and whether this procedure, without any theoretical approximation, is accurate enough. We have found that this approximation is very accurate for any value of the potential width and of the reduced temperature at high densities, except perhaps for small potential widths at very low temperatures. For lower densities the convergence slows down more markedly the lower are the values of the potential width and/or the reduced density. In addition, we have compared the theoretical predictions for the first-order contribution obtained from the integral equation theory of Tang and Lu (TL) [J. Chem. Phys. 100 (1994) 6665] with the simulation data. We have found that the theory is fairly accurate for distances larger than the width of the potential, whereas for shorter distances the accuracy decreases as the width of the potential decreases. (C) 2003 Elsevier B.V. All rights reserved.