  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.
