The ocean waves exhibit obvious non-linearity with asymmetric distribution of wave crests and troughs, which could induce significantly different effect on the seabed compared to the commonly used linear wave theory. In this paper, a semi-analytical solution for a transversely isotropic and multilayered poroelastic seabed under non-linear ocean wave is proposed by virtue of the dual variable and position (DVP) method. The ocean wave and seabed are, respectively, modelled using second-order Stokes theory and Biot’s complete poroelastodynamic theory. Then the established governing equations are decoupled and solved via the powerful scalar potential functions. Making use of the DVP scheme, the layered solutions are finally gained by combining the boundary conditions of the seabed. The developed solutions are verified by comparing with existing solutions. The selected numerical examples are presented to investigate the effect of main parameters on the dynamic response of the seabed and evaluate the corresponding liquefaction potential. The results show that the anisotropic stiffness and permeability, degree of saturation and stratification have remarkable influence on the dynamic response and liquefaction behavior of the seabed. The present solution is a useful tool to estimate the stability of transversely isotropic and layered seabed sediment in the range of non-linear ocean wave.