In this work, the planetary averaged equations are computed to high order using dedicated computer algebra. The resulting system comprises more than 150 000 polynomial terms in complex variables and approximate very well the long term evolution of the Solar System. The system of equations is developed up to order 2 in the masses and order 5 in eccentricity and inclination. The system of degree 3 is analytically integrated at first order and provide a solution with more than 25 000 terms. The problem of secular small divisors is discussed and a list of small divisors that arise with large amplitude in the solution is exhibited, the leading one being related to $g_1 -g_5 +(s_2 -s_1)$ that arises in the eccentricity of Mercury and Jupiter, and in the inclination of Mercury and Venus. It is acknowledge that the presence of these divisors greatly compromise the construction of the analytical solution, without considering the fact that the same solution of degree 5 comprises more than 3 000 000 terms. The seculars equations are then numerically integrated in a very efficient way over 1 million years, using a step size of 500 years for all the planets, after adding the averaged contribution from general relativity and the Moon. By comparison with the numerical ephemeris DE102 over more than 3000 years, it is demonstrated that the secular equations represent very well the long term evolution of the Solar System.
