This paper compares the efficiency and accuracy of different second-order approximations for the single-step propagator obtained by means of the power series expansion formalism and the split operator method. The former is based on a harmonic reference system, while the latter allows for the use of physically motivated (anharmonic) zeroth-order representations. Three typical examples---a system-bath Hamiltonian, a H\'enon-Heiles anharmonic resonating system, and a Fokker-Planck chaotic model---are considered in the present testing. The examples cover a variety of situations with strongly anharmonic coupling. Although no system-specific reference system is involved in the power series representation of the propagator, it quite accurately describes the dynamics of very anharmonic processes in the entire time domain even though the coupling is strong. Another appealing feature of the approach is that it is essentially analytical and therefore does not require any computational effort to implement. This makes the power series expansion technique particularly attractive for efficiently treating many-body problems. In contrast, numerical implementation of the split operator method can be arduous as a general multidimensional calculation, while its utility is in general restricted to the separable limit, when the coupling is almost turned off. In that limit the method indeed provides accurate results over a broad range of t. With increasing coupling, however, the efficiency of the method deteriorates very rapidly regardless of the particular choice of the zeroth-order representation.