  We present a method for computing actions of the exponential-like
$\varphi$-functions for a Kronecker sum $K$ of $d$ arbitrary matrices $A_\mu$.
It is based on the approximation of the integral representation of the
$\varphi$-functions by Gaussian quadrature formulas combined with a scaling and
squaring technique. The resulting algorithm, which we call PHIKS, evaluates the
required actions by means of $\mu$-mode products involving exponentials of the
small sized matrices $A_\mu$, without forming the large sized matrix $K$
itself. PHIKS, which profits from the highly efficient level 3 BLAS, is
designed to compute different $\varphi$-functions applied on the same vector or
a linear combination of actions of $\varphi$-functions applied on different
vectors. In addition, thanks to the underlying scaling and squaring techniques,
the desired quantities are available simultaneously at suitable time scales.
All these features allow the effective usage of PHIKS in the exponential
integration context. In fact, our newly designed method has been tested on
popular exponential Runge--Kutta integrators of stiff order from one to four,
in comparison with state-of-the-art algorithms for computing actions of
$\varphi$-functions. The numerical experiments with discretized semilinear
evolutionary 2D or 3D advection--diffusion--reaction, Allen--Cahn, and
Brusselator equations show the superiority of the proposed $\mu$-mode approach.
