We present the theoretical framework to efficiently solve the Jeans equations for multi-component axisymmetric stellar systems, focusing on the scaling of all quantities entering them. The models may include an arbitrary number of stellar distributions, a dark matter halo, and a central supermassive black hole; each stellar distribution is implicitly described by a two- or three-integral distribution function, and the stellar components can have different structural (density profile, flattening, mass, scale-length), dynamical (rotation, velocity dispersion anisotropy), and population (age, metallicity, initial mass function, mass-to-light ratio) properties. In order to determine the ordered rotational velocity and the azimuthal velocity dispersion fields of each component, we introduce a decomposition that can be used when the commonly adopted Satoh decomposition cannot be applied. The scheme developed is particularly suitable for a numerical implementation; we describe its realisation within our code JASMINE2, optimised to maximally exploit the scalings allowed by the Poisson and the Jeans equations, also in the post-processing procedures. As applications, we illustrate the building of three multi-component galaxy models with two distinct stellar populations, a central black hole, and a dark matter halo; we also study the solution of the Jeans equations for an exponential thick disc, and for its multi-component representation as the superposition of three Miyamoto-Nagai discs. A useful general formula for the numerical evaluation of the gravitational potential of factorised thick discs is finally given.