Exponential integrators require a reliable and efficient implementation of the action of the matrix exponential and related $\varphi$ functions. In this work we consider the Leja method for performing this task. We give an overview on recent developments with a closer look on a new a posteriori error estimate. A numerical comparison of a revised implementation to other methods from the literature is presented. For the new error estimate we define the notion of a residual based estimate where the residual is obtained from differential equations defining the $\varphi$ functions. For the numerical investigation of the newly defined error estimate as well as the comparison of the Leja interpolation we rely on test examples from spatial discretizations of time dependent partial differential equations. The experiments show that our new error estimate is robust and that the Leja interpolation performs very well in comparison to other methods.