In this paper we propose a numerical method for approximating the product of a matrix function with multiple vectors by Krylov subspace methods combined with a QR decomposition of these vectors. This problem arises in the implementation of exponential integrators for semilinear parabolic problems. We will derive reliable stopping criteria and we suggest variants using up- and downdating techniques. Moreover, we show how Ritz vectors can be included in order to speed up the computation even further. By a number of numerical examples, we will illustrate that the proposed method will reduce the total number of Krylov steps significantly compared to a standard implementation if the vectors correspond to the evaluation of a smooth function at certain quadrature points. Key words: Krylov subspace methods, shift and invert Lanczos process, projection method, matrix exponential function, matrix functions, restarts, error analysis, QR decomposition, multiple right-hand sides, exponential integ...