  Context. Dynamically self-consistent galactic models are necessary for
analysing and interpreting star counts, stellar density distributions, and
stellar kinematics in order to understand the formation and the evolution of
our Galaxy.
  Aims. We modify and improve the dynamical self-consistency of the
Besan\c{c}on Galaxy model in the case of a stationary and axisymmetric
gravitational potential.
  Methods. Each stellar orbit is modelled by determining a Staeckel approximate
integral of motion. Generalised Shu distribution functions with three integrals
of motion are used to model the stellar distribution functions.
  Results. This new version of the Besan\c{c}on model is compared with the
previous axisymmetric BGM2014 version and we find that the two versions have
similar densities for each stellar component. The dynamically self-consistency
is improved and can be tested by recovering the forces and the potential
through the Jeans equations applied to each stellar distribution function.
Forces are recovered with an accuracy better than one per cent over most of the
volume of the Galaxy.
