Although Kepler’s laws can be empirically proven by applying Newton’s laws to the dynamics of two particles attracted by gravitational interaction, an explicit formula for the motion as a function of time remains undefined. This paper proposes a quasi-analytical solution to address this challenge. It approximates the real dynamics of celestial bodies with a satisfactory degree of accuracy and minimal computational cost. This problem is closely related to Kepler’s equation, as solving the equations of motion as a function of time also provides a solution to Kepler’s equation. The results are presented for each planet of the solar system, including Pluto, and the solution is compared against real orbits.