We derive an algorithm for computing the wave-kernel functions and for an arbitrary square matrix , where . The algorithm is based on Padé approximation and the use of double angle formulas. We show that the backward error of any approximation to can be explicitly expressed in terms of a hypergeometric function. To bound the backward error we derive and exploit a new bound for that is sharper than one previously obtained by Al-Mohy and Higham [SIAM J. Matrix Anal. Appl., 31 (2009), pp. 970--989]. The amount of scaling and the degree of the Padé approximant are chosen to minimize the computational cost subject to achieving backward stability for in exact arithmetic. Numerical experiments show that the algorithm behaves in a forward stable manner in floating-point arithmetic and is superior in this respect to the general purpose Schur--Parlett algorithm applied to these functions.