A set of wave equations derived on the basis of a variational principle in consideration of both strong nonlinearity and strong dispersion of surface/internal waves is numerically solved to simulate generation and propagation of tsunamis in the vertical two-dimension. The velocity potential in each fluid layer is expanded into a power series of vertical position, such that the accuracy of vertical distribution of velocity depends on the number of expansion terms. Numerical results of surface displacement are compared with the existing experimental data, where tsunamis are generated by the seabed uplift. When the fundamental equations are reduced to nonlinear shallow water equations, the numerical model cannot represent propagation of a long wave group especially in distant-tsunami propagation, leading to overestimation of both the wave height and wave steepness of the first wave. The wave height becomes larger in the stratified ocean than that in a one-layer case, although the present density distribution hardly affects the tsunami phase over a long-distance travel.