  In this paper, we explore quadratures for the evaluation of $B^T \phi(A) B$
where $A$ is a symmetric nonnegative-definite matrix in $\mathbb{R}^{n \times
n}$, $B$ is a tall matrix in $\mathbb{R}^{n \times p}$, and $\phi(\cdot)$
represents a matrix function that is regular enough in the neighborhood of
$A$'s spectrum, e.g., a Stieltjes or exponential function. These formulations,
for example, commonly arise in the computation of multiple-input
multiple-output (MIMO) transfer functions for diffusion PDEs.
  We propose an approximation scheme for $B^T \phi(A) B$ leveraging the block
Lanczos algorithm and its equivalent representation through Stieltjes matrix
continued fractions. We extend the notion of Gauss-Radau quadrature to the
block case, facilitating the derivation of easily computable error bounds. For
problems stemming from the discretization of self-adjoint operators with a
continuous spectrum, we obtain sharp estimates grounded in potential theory for
Pad\'e approximations and justify extrapolation algorithms at no added
computational cost. The obtained results are illustrated on large-scale
examples of 2D diffusion and 3D Maxwell's equations as well as a graph from the
SNAP repository. We also present promising experimental results on convergence
acceleration via random enrichment of the initial block $B$.
