We propose a new approach for computing the action of the matrix exponential over a vector. The approach is based on a Legendre polynomial expansion in the framework of the so‐called ‐product solution to ODEs. The new approach can be combined with Krylov subspace methods.