Several authors have proposed studying randomly forced turbulent flows (e.g., E. A. Novikov, Soviet Physics JETP, 20(5), 1290 1965). More recently, theoretical investigations (e.g., renormalization group) have focused on white-noise forced Navier-Stokes equations (V. Yakhot and S. A. Orszag, J. Sci. Comput. 1(1), 3 1986). The present article aims to provide an appropriate numerical method for the simulation of randomly forced turbulent systems. The spatial discretization is based on the classical Fourier spectral method. The time integration is performed by a secondorder Runge-Kutta scheme. The consistency of an extension of this scheme to treat additive noise stochastic differential equations is proved. The random number generator is based on lagged Fibonacci series. Results are presented for two randomly forced problems: the Burgers and the incompressible Navier-Stokes equations with a white-noise in time forcing term characterized by a power-law correlation function in spectral space. A variety of statistics are computed for both problems, including the structure functions. The third-order structure functions are compared with their exact expressions in the inertial subrange. The influence of the dissipation mechanism (viscous or hyperviscous) on the inertial subrange is discussed. In particular, probability density functions of velocity increments are computed for the Navier-Stokes simulation. Finally, for both Burgers and Navier-Stokes problems, our results support the view that random sweeping is the dominant effect of the large-scale motion on the small-scales.