In linear, three-dimensional, continuously stratified equatorial β-plane ocean models with arbitrary eastern and western boundaries the shallow water equations for each vertical mode must be solved numerically in the horizontal variables. This paper introduces a new numerical method of solution for the time-Fourier transformed shallow water equations with slip boundary conditions at boundaries of arbitrary geometry. The method is based on a boundary integral equation (BIE) for the pressure perturbation response to the specified wind-stress forcing field. All other dependent variables are expressed as boundary functionals of the pressure perturbation. The kernels of all functionals are constructed from the Green's function for the Laplace Tidal Equation on the β-plane and its derivatives. The efficient computation of these kernels from their exact meridional mode representations may be performed by use of asymptotic methods especially developed for the numerical evaluation of functions expressed as slowly converging series of Hermite functions. The solution of the basic BIE and the computation of the boundary functionals involve the discretization of the ocean boundaries into a number of boundary segments (boundary elements). It is shown that the terms of the BIE involving the wind-stress field over the ocean may be reduced to a boundary integral, which effectively reduces the simulation problem to the solution of a one-dimensional BIE. The method incorporates the ocean physics through the relationship between the coastal pressure field and the basin-wide variables, pointing out to the possibility that the dynamic topography of the ocean may be estimated directly from the wind-stress field and coastal sea-level data.