In this study a general formula was established for multi-center integrals over Gaussian basis sets using the two-center expansion of nuclear attraction and electron repulsion potentials. The terms of the formula are the analytical functions of interatomic distances. Three kinds of integrals contained in these formulae can be evaluated analytically. In the presented method, unlike many other popular methods, the numerical integrations procedure is unnecessary since it does not contain the incomplete gamma function. A computer program was written and tested for various orbitals. The calculations show that the series converge very rapidly.