We have studied the dispersion relation for surface plasmons on a randomly rough surface, going beyond the lowest approximation in the surface profile function. The expansion of the surface plasmon dispersion relation in powers of the surface profile function contains an infinite subset of terms that are all of the same order of magnitude as the lowest order contribution, the only one considered in previous theoretical determinations of this relation. This subset of terms yields a nonlinear integral equation for the surface plasmon proper self-energy, in terms of which the dispersion relation is expressed. The splitting of the surface plasmon dispersion curve into two branches by the surface roughness, predicted theoretically by the lowest order perturbation theory calculation, and observed experimentally, is preserved in the results of the present calculation. However, both the magnitude of the splitting and the damping of the surface plasmon obtained here are larger, for the same corrugation strength, than the same quantities obtained from the lowest order perturbation calculation, for most surface plasmon wave vectors.