Abstract High-precision ephemerides not only support space missions, but can also be used to study the origin and future of celestial bodies. In this paper, a coupled orbit-rotation dynamics model that fully takes into account the rotation of the Martian moons is developed. Phobos and Deimos' rotation are firstly described by Eulerian rotational equations, and integrated simultaneously with the orbital motion equations. Orbital and orientational parameters of Mars satellites were simultaneously obtained by numerical integration for the first time. In order to compare the differences between our newly developed model and the one now used in the ephemerides, we first reproduced and simulated the current model using our own parameters, and then fit it to the IMCCE ephemerides using least-square procedures. The adjustment test simulations show Phobos and Deimos‘ orbital differences between the refined model and the current model is no more than 300~$m$ and 125~$m$, respectively. The orientation parameters are confirmed and the results are in good agreement with the IAU results. Moreover, we simulated two perturbations (main asteroids and mutual torques) which were not included in our refined model, and find that their effects on the orbits are completely negligible. As for the effect on rotation, we propose to take care of the role of mutual attraction in future models.