We develop a robust numerical method for modelling nonlinear gravity waves which is based on the Zakharov equation/mode-coupling idea but is generalized to include interactions up to an arbitrary order M in wave steepness. A large number (N = O(1000)) of free wave modes are typically used whose amplitude evolutions are determined through a pseudospectral treatment of the nonlinear free-surface conditions. The computational effort is directly proportional to N and M, and the convergence with N and M is exponentially fast for waves up to approximately 80% of Stokes limiting steepness (ka ∼ 0.35). The efficiency and accuracy of the method is demonstrated by comparisons to fully nonlinear semi-Lagrangian computations (Vinje & Brevig 1981); calculations of long-time evolution of wavetrains using the modified (fourth-order) Zakharov equations (Stiassnie & Shemer 1987); and experimental measurements of a travelling wave packet (Su 1982). As a final example of the usefulness of the method, we consider the nonlinear interactions between two colliding wave envelopes of different carrier frequencies.