  This report is devoted to the analytical calculation of heavy quarkonia
production processes in modern experiments such as LHC, B-factories and
superB-factories in computer. Theoretical description of heavy quarkonia is
based on the factorization theorem. This theorem leads to special structure of
the production amplitudes which can be used to develop computer algorithm which
calculates these amplitudes automatically. This report is devoted to the
description of this algorithm. As an example of its application we present the
results of the calculation of double charmonia production in bottomonia decays
and inclusive the $\chi_{cJ}$ mesons production in pp-collisions.
