A variety of block Krylov subspace methods have been successfully developed for linear systems and matrix equations. The application of block Krylov methods to computing matrix functions is, however, less established, despite the growing prevalence of matrix functions in scientific computing. Of particular importance is the evaluation of a matrix function on not just one but multiple vectors. The main contribution of this work is a class of efficient block Krylov subspace methods tailored precisely to this task. With the full orthogonalization method (FOM) for linear systems forming the backbone of our theory, the resulting methods are referred to as B(FOM)$^2$: block FOM for functions of matrices. Many other important results are obtained in the process of developing these new methods. Matrix-valued inner products are used to construct a general framework for block Krylov subspaces that encompasses already established results in the literature. Convergence bounds for B(FOM)$^2$... more