In this report, we study the system (*) N-t + W-x + (NW)(x) - 1/6N(xxt) = 0, W-t + N-x + WWx - 1/6W(xxt) = 0, which describes approximately the two-dimensional propagation of surface waves in a uniform horizontal channel of length L-0 filled with an irrotational, incompressible, inviscid fluid which in its undisturbed state has depth h. The non-dimensional variables N(x, t) and W(x, t) represent the deviation of the water surface from its undisturbed position and the horizontal velocity at water level root 2/3h, respectively. The natural initial-boundary-value problem corresponding to the situation wherein the channel is fitted with a wavemaker at both ends is formulated and analyzed theoretically, We then present a numerical algorithm for the approximation of solutions of the system (*) and prove the algorithm is fourth-order accurate both in time and in space, is unconditionally stable, and has optimal computational complexity, which is to say the operation cost per time step is of order M where M is the number of grid points in the spatial discretization. Further, we implement the algorithm as a computer code and use it to study head-on collisions of solitary waves. Our numerical simulations are compared with existing theoretical, numerical and experimental results. We have tentatively concluded that the system (*) is a good candidate for modeling two-dimensional surface water waves. Copyright (C) 1998 Elsevier Science B.V.