For the time evolution of acoustic wavefields we present an alternative derivation of the pseudo-analytical method, which enables us to generalize the method to high-order formulations. Within the same derivation framework, we compare the second-order pseudo-analytical method, the Fourier finite difference method, and the fourth-order Lax-Wendroff time integration method. We demonstrate that the pseudo-analytical method can be regarded as a modified Lax-Wendroff method. Different from the fourth-order time stepping method, both the second-order pseudo-analytical method and the Fourier finite difference method use pseudo-Laplacians to compensate for time stepping errors. The pseudo-Laplacians need to be solved in the wavenumber domain with constant compensation velocities for computational simplicity and efficiency. Low-order pseudo-Laplacians are more sensitive to the choice of compensation velocities than high-order ones. As a result, we need to use the combination of several pseudo-Laplacians to achieve the required accuracy for low-order pseudo-analytical methods. When using the pseudospectral method to evaluate all spatial derivatives, the computation cost for the second-order pseudo-analytical method, the Fourier finite difference method, and the fourth-order Lax-Wendroff time integration method is approximately the same. Both the second-order pseudo-analytical method and the Fourier finite difference method have less restrictive stability conditions than the fourth-order time stepping method. We demonstrate with numerical examples that the second-order pseudo-analytical method, greatly improves the original pseudo-analytical method and as a modified version of the Lax-Wendroff method, is well suited for imaging seismic data in subsalt areas where reverse-time migration plays a crucial role.