A spectral method is developed to numerically solve the so-called Kuramoto-Sakaguchi equation, which is a nonlinear integro-differential equation of the parabolic type, governing the dynamical statistical behaviour of certain populations of nonlinearly coupled random oscillators. The approach rests on explicit bounds for the space derivatives of solutions, obtained via energy-like estimates. Bounds for the numerical approximations of solutions are given, and improved (sometimes appreciably) by means of an ‘a posteriori error analysis’. Plots are shown to illustrate the performance of the method, and comparison with a finite difference approach is also made.