This paper presents a numerical approach based on spectral methods for the computation of guided ultrasonic wave modes (i.e., Lamb and shear horizontal) in nonuniformly stressed plates. In particular, anisotropic elastic plates subjected to a normal stress profile, which varies nonuniformly over their thickness, are considered. The proposed approach computes the modeshapes and the full three-dimensional dispersion spectrum (i.e., real frequency, complex wavenumber). It therefore includes both propagating (real wavenumber) and non-propagating (complex wavenumber) modes. Furthermore, an approach for robustly post-processing the dispersion spectra in order to compute the group velocity of propagating modes is presented, which is based on a spectral quadrature method. Numerical results are presented for two case studies: (1) a bending profile in a fiber-reinforced graphite/epoxy plate, and (2) an exponential profile in a silver plate. The results show the computational efficiency (i.e., spectral convergence) of the proposed method compared to other existing approaches such as the sublayering and finite element methods.