We consider a new Krylov subspace algorithm for computing expressions of the form $\sum\limits_{k=0}^p h^k \varphi_k(hA)w_k$, where $A\in \mathbb{R}^{n \times n}$, $w_k \in \mathbb{R}^n$, and $\varphi_k$ are matrix functions related to the exponential function. Computational problems of this form appear when applying exponential integrators to large dimensional ODEs in semilinear form $u'(t) = Au(t) + g(u(t))$. Using Cauchy's integral formula we give a representation for the error of the approximation and derive a priori error bounds which describe well the convergence behaviour of the algorithm. In addition an efficient a posteriori estimate is derived. Numerical experiments in MATLAB illustrating the convergence behaviour are given.