We present an efficient method for simulating a stationary Gaussian noise with an arbitrary covariance function, and then we study numerically the impact of time-correlated noise on the time evolution of a (1+1)-dimensional generalized Langevin equation by comparing also to analytical results. Finally, we apply our method to the generalized Langevin equation with an external harmonic and double-well potential.