The two-center overlap integrals in the molecular coordinate system over Slater-type orbitals are expressed in terms of newly defined functions, which are evaluated recursively and are closely related to the well-known incomplete gamma functions. Using the well-known rotational properties of spherical harmonics, the rotation matrices, also referred to as rotation coefficients, were calculated. Rotation matrices were expressed in terms of Gaunt coefficients and complex spherical harmonics. They were computed in two different ways: in the first approach, the Gaunt coefficients and normalized complex spherical harmonics were directly evaluated using binomial coefficients, while in the second, both quantities were obtained recursively. The numerical stability of our computational results was analyzed. In the calculation of overlap integrals and rotation matrices, both methods were compared in terms of accuracy and computation time (CPU) for low and very high quantum numbers.