HELAC-Onia: an automatic matrix element generator for heavy quarkonium physics. By the virtues of the Dyson-Schwinger equations, we upgrade the published code HELAC to be capable to calculate the heavy quarkonium helicity amplitudes in the framework of NRQCD factorization, which we dub HELAC-Onia. We rewrote the original HELAC to make the new program be able to calculate helicity amplitudes of multi P-wave quarkonium states production at hadron colliders and electron-positron colliders by including new P-wave off-shell currents. Therefore, besides the high efficiencies in computation of multi-leg processes within the Standard Model, HELAC-Onia is also sufficiently numerical stable in dealing with P-wave quarkonia (e.g. h c,b ,χ c,b ) and P-wave color-octet intermediate states. To the best of our knowledge, it is a first general-purpose automatic quarkonium matrix elements generator based on recursion relations on the market.