Over recent decades, fiber Raman lasers (FRLs) have received much attention from researchers and have become a challenge for them both numerically and experimentally. The equations governing the FRLs are in the form of a first-order system of nonlinear two-point boundary-value ordinary differential equations. In this paper, an algorithm for solving this system of differential equations using a spectral method, namely Chebyshev pseudospectral method, is presented in detail and then numerical simulations are performed. The main advantage of the spectral methods is in their optimality in achieving high accuracy by using fewer degrees of freedom under suitable conditions. It is shown that the proposed spectral method in combination with the Newton method results in a considerable reduction in the size of the discretized problem and in the computational effort to achieve high accuracy. In this paper, a new approach for constructing an initial approximate solution for the Newton iteration is also presented.