Non-hydrostatic multi-layer models have become a popular tool in describing wave transformation from deep water to the surf zone, but the numerical approach lacks a theoretical framework to guide implementation and assist interpretation of the results. In this paper, we formulate a non-hydrostatic model in an analytical form for the derivation and examination of dispersive and nonlinear properties. Depth integration of the dimensionless continuity and Euler equations over each layer yields the conventional multi-layer formulation. A variable transformation converts the conventional form into an integrated series form, which provides separate descriptions of flux- and dispersion-dominated processes. Substitution of the non-hydrostatic pressure and vertical velocity in the governing equations by high-order derivatives of the horizontal velocity and surface elevation provides a direct comparison with the Boussinesq equations published in the literature. Implementation of a perturbation expansion extracts the first- and second-order governing equations with respect to the nonlinear parameter. Based on that, we derive analytical solutions of the linear dispersion and the second-order super- and sub-harmonics for up to three layers and optimize the solutions in terms of the layer arrangement. In relation to the Boussinesq equations at comparable orders of expansion, the two- and three-layer models provide slightly higher errors in shallow and intermediate water in terms of dispersion and super-harmonics, but show superior performance in describing sub-harmonics in deep water.