Starting with an exact and simple geodesic, we generate approximate geodesics1,2 by summing up higher-order geodesic deviations in a fully relativistic scheme. We apply this method to the problem of orbit motion of test particles in Schwarzschild3 and Kerr metrics; from a simple circular orbit as the initial geodesic we obtain finite eccentricity orbits as a Taylor series with respect to the eccentricity. The explicit expressions of these higher-order geodesic deviations are derived using successive systems of linear equations with constant coefficients, whose solutions are of harmonic oscillator type. This scheme is best adapted for small eccentricities, but arbitrary values of M/R. We also analyse the possible application to the calculation of the emission of gravitational radiation from non-circular orbits around a very massive body3.