In the present work, a stable algorithm for the calculation of the electromagnetic field distributions of the eigenmodes of one-dimensional diffraction gratings is presented. The proposed approach is based on the method for the computation of the propagation constants of Bloch waves of such structures previously presented by Cao et al.[J. Opt. Soc. Am. A 19, 335 (2002)] and uses a modified S-matrix algorithm to ensure numerical stability.