An alternative methodology to evaluate two-electron-repulsion integrals based on numerical approximation is proposed. Computational chemistry has branched into two major fields with methodologies based on quantum mechanics and classical force fields. However, there are significant shadowy areas not covered by any of the available methods. Many relevant systems are often too big for traditional quantum chemical methods while being chemically too complex for classical force fields. Examples include systems in nanomedicine, studies of metalloproteins, etc. There is an urgent need to develop fast quantum chemical methods able to study large and complex systems. This work is a proof-of-concept on the numerical techniques required to develop accurate and computationally efficient algorithms for the fast calculation of electron-repulsion integrals, one of the most significant bottlenecks in the extension of quantum chemistry to large systems. All concepts and calculations were developed for the three-center integral (p xA p xB |p xC p xC ) with all atoms being carbon. Starting with the analytical formulae, convenient decompositions were tested to provide smooth two-dimensional surfaces that were easily fitted. The approximating algorithm consisted of a multilayered approach based on multiple fittings of two-dimensional surfaces. An important aspect of the new method is its independence on the number of contracted Gaussian primitives. The basis set of choice was STO-6G. In future developments, larger basis sets will be developed. This work is part of a large effort aimed at improving the inadequacies of existing computational chemistry methods, both based on quantum mechanics and classical force fields, in particular in describing large and heterogeneous systems (ex. metalloproteins).