A new method valid for highly dispersive and highly nonlinear water waves is presented. It combines a time-stepping of the exact surface boundary conditions with an approximate series expansion solution to the Laplace equation in the interior domain. The starting point is an exact solution to the Laplace equation given in terms of infinite series expansions from an arbitrary z-level. We replace the infinite series operators by finite series (Boussinesq-type) approximations involving up to fifth-derivative operators. The finite series are manipulated to incorporate Padé approximants providing the highest possible accuracy for a given number of terms. As a result, linear and nonlinear wave characteristics become very accurate up to wavenumbers as high as kh = 40, while the vertical variation of the velocity field becomes applicable for kh up to 12. These results represent a major improvement over existing Boussinesq-type formulations in the literature. A numerical model is developed in a single horizontal dimension and it is used to study phenomena such as solitary waves and their impact on vertical walls, modulational instability in deep water involving recurrence or frequency downshift, and shoaling of regular waves up to breaking in shallow water.