Context. Over short time-intervals, planetary ephemerides have traditionally been represented in analytical form as finite sums of periodic terms or sums of Poisson terms that are periodic terms with polynomial amplitudes. This representation is not well adapted for the evolution of planetary orbits in the solar system over million of years which present drifts in their main frequencies as a result of the chaotic nature of their dynamics. Aims. We aim to develop a numerical algorithm for slowly diffusing solutions of a perturbed integrable Hamiltonian system that will apply for the representation of chaotic planetary motions with varying frequencies. Methods. By simple analytical considerations, we first argue that it is possible to exactly recover a single varying frequency. Then, a function basis involving time-dependent fundamental frequencies is formulated in a semi-analytical way. Finally, starting from a numerical solution, a recursive algorithm is used to numerically decompose the solution into the significant elements of the function basis. Results. Simple examples show that this algorithm can be used to give compact representations of different types of slowly diffusing solutions. As a test example, we show that this algorithm can be successfully applied to obtain a very compact approximation of the La2004 solution of the orbital motion of the Earth over 40 Myr ([−35 Myr, 5 Myr]). This example was chosen because this solution is widely used in the reconstruction of the past climates.