Computing Linear Combinations of $\varphi$-Function Actions for Exponential Integrators

We propose a matrix-free algorithm for evaluating linear combinations of $\varphi$-function actions, $w_i := \sum_{j=0}^{p} \alpha_i^{\,j}\,\varphi_j(t_i A)v_j$ for $i=1\colon r$, arising in exponential integrators. The method combines the scaling and recovering method with a truncated Taylor series, choosing a spectral shift and a scaling parameter by minimizing a power-based objective of the shifted operator. Accuracy is user-controlled and ultimately limited by the working precision. The algorithm decouples the stage abscissae $t_i$ from the polynomial weights $\alpha_i^j$, and a block variant enables simultaneous evaluation of $\{w_i\}_{i=1}^r$. Across standard benchmarks, including stiff and highly nonnormal matrices, the algorithm attains near-machine accuracy (IEEE double precision in our tests) for small step sizes and maintains reliable accuracy for larger steps where several existing Krylov-based algorithms deteriorate, providing a favorable balance of reliability and computational cost.