The idea is to use holomorphic functional calculus and to show that for a matrix $A$ and a polynomial $p(z)$ we have for $r \gt \|A\|$
\begin{equation}\tag{$\ast$}
p(A) = \frac{1}{2\pi i} \int_{|z| = r} p(z) \cdot (z - A)^{-1}\ \,dz
\end{equation}
in complete analogy with the Cauchy formula for complex numbers. The integral of a matrix of holomorphic functions is defined by integrating each entry separately.
By Cramer's rule, the $(k,l)$-entry of $(z-A)^{-1}$ is $\displaystyle ((z-A)^{-1})_{k,l} = \frac{1}{\det(z-A)} c_{k,l}(z)$ where $c_{k,l}(z)$ is some polynomial in $z$. Let $p(z) = \det(z-A)$ be the characteristic polynomial of $A$. Conclude using $(\ast)$ by applying Cauchy's integral theorem to $c_{k,l}$.
To see that the identity $(\ast)$ holds, proceed as follows (this is a slight variant of McCarthy's argument):
- The usual matrix norm induced by the Euclidean norm on $\mathbb{C}^{n}$ satisfies $\|A^{n}\| \leq \|A\|^{n}$.
- Use this to show that $(z - A)^{-1} = \sum_{n = 0}^{\infty} \frac{A^{n}}{z^{n+1}}$, where the right hand side converges uniformly on $\{|z| \gt \|A\| + \varepsilon\}$.
- It follows that we can interchange integration and summation. Conclude that $$ A^{k} = \int_{|z| = r} z^{k} (z - A)^{-1}\,dz$$ and $(\ast)$ follows by linearity.
Here's a link to McCarthy's article (you need a university subscription to download it, but the first page is almost the entire article).