We present an original method for integrating numerically a scalar generalized Langevin Equation (GLE) for an arbitrary gaussian ‘random force’ process. Our procedure rests on the approximation of this last stochastic process by a vectorial gaussian-Markov one which reproduces the original autocorrelation function within a given accuracy. This is obtained by truncating at some suitablenth order the Mori continued fraction representation of the original autocorrelation function. In this way, the integrodifferential GLE is shown to be equivalent to a set ofnfirst order differential equations with an additional gaussian white noise source. This system can then be integrated with the usual algorithms. The method is illustrated by the generalized brownian motion of a spherical particle in a medium at equilibrium. Assuming a gaussian velocity autocorrelation function, we compare Markov approximations of increasing order to the original function and we perform a stochastic dynamics of the process in the 10th order Markov approximation. Various time averages computed on the resulting trajectory are shown to be in good agreement with the known exact results.