A new analytic solution is presented for tidal propagation in multilayer coastal aquifers consisting of an arbitrary number of layers. The solution is derived using matrix calculus and may be applied to simulate tidal propagation in systems where aquifer layers are separated by leaky layers of low permeability or in stratified aquifers. The head, amplitude, and phase shift may be simulated for arbitrary tidal components. The transmissivity, storage coefficient, and loading efficiency may be different for every aquifer layer. Similarly, the resistance to vertical flow, storage coefficient, and loading efficiency may be different for every leaky layer. The aquifer system extends an infinite distance below the sea or may end abruptly at the shore line. In the former case, different aquifer and leaky layer properties may be specified below the sea and below the land. A homogeneous unconfined aquifer is simulated with 80 layers to demonstrate that the tidal signal propagates much farther inland in the deeper part than in the shallower part of an unconfined aquifer. This effect is enhanced significantly when layers of low permeability are present in the aquifer, even when these layers are fairly thin, such as clay lenses. The solution is implemented in Python and may be used, for example, to estimate aquifer parameters from head measurements in coastal aquifers.