A numerical method is developed to solve a class of nonlinear, nonlocal eigenvalue problems defined in an infinite strip, and is applied to compute solitary planetary waves in a sheared zonal current on the beta‐plane. This method, an iterative procedure derived from the natural variational structure of these problems, is implemented in the physical case when the ambient parallel flow has a linear or a quadratic velocity profile. The results of the numerical experiments establish rigorous limits on the range of validity of the formal asymptotic theory of weakly nonlinear long waves, and also reveal some new phenomena involving strongly nonlinear waves. The iterative procedure is analyzed in a general setting, and is shown to be globally convergent without restriction on the wave amplitud