This article is concerned with the construction of the general algorithm for evaluating twocenter, two-and three-electron integrals occurring in matrix elements of one-electron operators in the basis of variational correlated functions. This problem has been solved here in prolate spherical coordinates, using the modified and extended form of the Neumann expansion of the interelectronic distance function rt derived in Part I of this series for k = -1,O, 1,2. This work expands the method proposed by one of us in the preceding paper for integrals of the types mentioned above. The results of numerical calculaticns for different types of the two-and threeelectron integrals are presented. The problem of convergence of the proposed procedures used is also discussed.