2026-07-072026-07-07http://salesiana.dossiersoluciones.com/handle/123456789/100707A method for the evaluation of the epsilon expansion of multi-loop massless Feynman integrals is introduced. This method is based on the Gegenbauer polynomial technique and the expansion of the Gamma function in terms of harmonic sums. Algorithms for the evaluation of nested and harmonic sums are used to reduce the expressions to get analytical or numerical results for the expansion coefficients. Methods to increase the precision of numerical results are discussed.30 pages, 6 figures; Minor typos corrected, references added. Published in Computer Physics CommunicationsHigh Energy Physics - Phenomenologytext