Let $A$ be a square complex matrix; $z_1$, ..., $z_{N}inmathbb C$ be arbitrary (possibly repetitive) points of interpolation; $f$ be an analytic function defined on a neighborhood of the convex hull of the union of the spectrum $sigma(A)$ of the matrix $A$ and the points $z_1$, ..., $z_{N}$; and the rational function $r=frac uv$ (with the degree of the numerator $u$ less than $N$) interpolates $f$ at these points (counted according to their multiplicities). Under these assumptions estimates of the kind $$ biglVert f(A)-r(A)bigrVertle max_{tin[0,1];muintext{convex hull}{z_1,z_{2},dots,z_{N}}}bigglVertOmega(A)[v(A)]^{-1} frac{bigl(vfbigr)^{{(N)}} bigl((1-t)mumathbf1+tAbigr)}{N!}biggrVert, $$ where $Omega(z)=prod_{k=1}^N(z-z_k)$, are proposed. As an example illustrating the accuracy of such estimates, an approximation of the impulse response of a dynamic system obtained using the reduced-order Arnoldi method is considered, the actual accuracy of the approximation is compared with the estimate based on this paper.