Over short time intervals planetary ephemerides have been traditionally represented in analytical form as finite sums of periodic terms or sums of Poisson terms that are periodic terms with polynomial amplitudes. Nevertheless, this representation is not well adapted for the evolution of the planetary orbits in the solar system over million of years as they present drifts in their main frequencies, due to the chaotic nature of their dynamics. The aim of the present paper is to develop a numerical algorithm for slowly diffusing solutions of a perturbed integrable Hamiltonian system that will apply to the representation of the chaotic planetary motions with varying frequencies. By simple analytical considerations, we first argue that it is possible to recover exactly 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 on the significant elements of the function basis. 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 how 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 ([-35Myr,5Myr]). This example has been chosen as this solution is widely used for the reconstruction of the climates of the past.