By utilizing a novel quasi-harmonic three-layer dielectric model for the interface between a dielectric spheroid and the surrounding dissimilar dielectric medium, a robust numerical method for calculating the generalized Coulomb and self-polarization potentials of the dielectric spheroid is presented in this paper. The proposed numerical method can not only overcome the inherent mathematical divergence in the self-polarization energy which arises for the simplest step-like model of the dielectric interface, but also completely eliminate the potential numerical divergence which may occur in other treatments. Numerical experiments have demonstrated the convergence of the proposed numerical method as the number of the steps used to discretize the translation layer in a general three-layer dielectric model goes to infinity.