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.