We investigate the use of a restarted Krylov subspace method for computing the action of a function of large sparse matrix on a vector to numerically solve the one-step wave extrapolation matrix (OSEM) method. The OSEM method belongs to a category exempt from stability and dispersion problems when extrapolating acoustic wavefields in time. It is based on the solution of a first-order differential equation in the time domain for the analytical wavefield. Such wavefield is separated into real and imaginary parts, and the resulting first-order matrix differential equation is integrated in time using Chebyshev polynomial expansion, usually combined with the Fourier method to compute a square-root pseudodifferential operator. To by- pass fast fourier transforms (FFT), the aim of the present work is to introduce an approach that uses the theory of matrix functions, and therefore enables the use of fast Krylov-based methods for their computation. Numerical modeling examples are shown to demonstrate the efficacy of the new formulation.