This paper shows how a self-consistent dynamical model can be obtained by fitting the gravitational potential of the Milky Way to the stellar kinematics and densities from Gaia data. Using the Besancon Galaxy Model we derive a potential and the disc stellar distribution functions are computed based on three integrals of motion to model stationary stellar discs. The gravitational potential and the stellar distribution functions are built self-consistently, and then adjusted to be in agreement with the kinematics and the density distributions obtained from Gaia observations. A Markov chain Monte Carlo (MCMC) is used to fit the free parameters of the dynamical model to Gaia parallax and proper motion distributions. The fit is done on several sets of Gaia eDR3 data, widely spread in longitudes and latitudes. We are able to determine the velocity dispersion ellipsoid and its tilt for sub-components of different ages, both varying with R and z. The density laws and their radial scale lengths, for the thin and thick disc populations are also obtained self-consistently. This new model has some interesting characteristics, such as a flaring thin disc. The thick disc is found to present very distinctive characteristics from the old thin disc, both in density and kinematics. This well supports the idea that thin and thick discs were formed in distinct scenarios as the density and kinematics transition between them is found to be abrupt. The dark matter halo is shown to be nearly spherical. We also derive the Solar motion to be (10.79 $\pm$ 0.56, 11.06 $\pm$ 0.94, 7.66 $\pm$ 0.43) km/s, in good agreement with recent studies. The resulting fully self-consistent gravitational potential, still axisymmetric, is a good approximation of a smooth mass distribution in the Milky Way and can be used for further studies, including to compute orbits for real stars in our Galaxy (abridged).