Solving integral equations is an effective approach to obtain the radial distribution function (RDF) of multicomponent mixtures. In this work, by extending Gillan's approach [M. J. Gillan, Mol. Phys. 38(6), 1781 (1979)], the integral equation was solved by numerical method and was applied to both binary and ternary mixtures. The Lennard-Jones (LJ) potential function was used to express the pair molecular interactions in calculating the RDF and chemical potential. This allowed a comparison with available simulation data, on the RDF and the chemical potential, since the simulation data have been reported for the LJ potential function. The RDF and the chemical potential results indicated good agreement with the simulation data. The calculations were extended to the ternary system and the RDFs for carbon dioxide-octane-naphthalene were obtained. The numerical method used in solving integral equation was rapidly convergent and not sensitive to the first estimation. The method proposed in this work can be easily extended to more than the three-component systems.