For a symmetric matrix $A$ with eigenvalues $\lambda_i<1$ (with $i=0,\dots,n$), one can write the resolvent ($R(z,A):=(A-zI)^{-1}$) at $z=1$ as a series of powers of $A$ $[1]$:
$R(1,A)=(A-1)^{-1}=-(I+A+A^2+A^3+\dots)$
From Cayley-Hamilton, I know that one can express the power $k\geq n$ of an $n\times n$ matrix $A$ as a linear combination of $A^0,A^1,A^2,\dots, A^{n-1}$. I would like to use this to convert the infinite series of matrix powers above into a finite series of matrix powers (with possibly very ugly coefficients). For this, I need a formula for the $k$-th power of an $n\times n$ matrix of the form
$A^k=\sum_{j=0}^{n-1}\xi_jA^j$
for $k\geq n$. The $\xi_j=\xi_j(c_0,c_1,c_2,\dots,c_n)$ should be function s of the coefficients $c_i$ (with $i=0,\dots,n$) of the $A$'s characteristic polynomial. I have tried to derive this formula by hand and so far it seems to be possible, but very tedious. I would assume that other people have worked on finding such a formula before, but I couldn't find it on the internet.
$[1]$ Barnett, Lionel, Christopher L. Buckley, and Seth Bullock. "Neural complexity and structural connectivity." Physical Review E 79.5 (2009): 051914.
Obviously the expansion must be unique, so we can find it however we like. If $$p(t) = \sum_{k=0}^n c_k t^k $$ is the characteristic polynomial of $A$, then if we write $B=A-I$, then $A=I+B$ satisfies $$ \sum_{k=0}^n c_k (I+B)^k = 0 $$ by Cayley–Hamilton. We can then expand to get $$ 0 = \sum_{k=0}^n c_k \sum_{m=0}^k \binom{k}{m} B^m = \sum_{m=0}^n B^m \sum_{k=m}^n c_k \binom{k}{m} = \sum_{m=0}^n \alpha_m B^m $$ by changing the order of summation, and the $\alpha_m = \sum_{k=m}^n c_k \binom{k}{m} $ are constants. So $$ \alpha_0 = -\sum_{m=1}^n B^m \alpha_m. $$ $\alpha_0$ is nonzero since $B$ is supposed to be invertible, so multiplying both sides by $B^{-1}\alpha^{-1}$, we have $$ B^{-1} = \sum_{m=1}^n -\frac{\alpha_m}{\alpha_0} B^{m-1} . $$ Now we can go back to $A$ using $B=A-I$, which gives $$ (A-I)^{-1} = \sum_{m=0}^{n-1} -\frac{\alpha_{m+1}}{\alpha_0} (A-I)^{m} = \sum_{m=0}^{n-1} -\frac{\alpha_{m+1}}{\alpha_0} \sum_{j=0}^m (-1)^{m-j} \binom{m}{j} A^j \\ = \sum_{j=0}^{n-1} A^j \frac{-1}{\alpha_0}\sum_{m=j}^{n-1} (-1)^{m-j} \binom{m}{j} \alpha_{m+1} \\ = \sum_{j=0}^{n-1} A^j \frac{-1}{\alpha_0}\sum_{m=j}^{n-1} (-1)^{m-j} \binom{m}{j} \sum_{k=m+1}^n c_k \binom{k}{m+1}. $$ Therefore $$ (A-I)^{-1} = \sum_{j=0}^{n-1} \beta_j A^j, $$ with coefficients $$ \beta_j = -\left( \sum_{k=0}^n c_k \right)^{-1} \sum_{m=j}^{n-1} (-1)^{m-j} \binom{m}{j} \sum_{k=m+1}^n c_k \binom{k}{m+1} $$ It may be possible to simplify this expression further by changing the order of summation again and using some binomial identities.
Edit:
$$ \sum_{m=j}^{n-1} (-1)^{m-j+1} \binom{m}{j} \sum_{k=m+1}^n c_k \binom{k}{m+1} = (-1)^j \sum_{k=j+1}^n c_k \sum_{M=j+1}^{k} (-1)^{M} \binom{M-1}{j} \binom{k}{M} = -\sum_{k=j+1}^n c_k $$ by Euler's formula on sums of $k$th differences of powers, since $j<k$.
Actually, since everything commutes, we can do some formal matrix division to get this result more simply: $$ 0 = (A-I)^{-1}(c_n A^n+ c_{n-1} A^{n-1} + \dotsb + c_1 A + c_0 I ). $$ Then $$ (A-I)^{-1}A^k = A^{k-1}+A^k+\dotsb+I + (A-I)^{-1}, $$ and so we have $$ \left( \sum_{k=0}^n c_k \right) (A-I)^{-1} = -\sum_{m=0}^{n} A^{m} \sum_{k=m+1}^n c_k, $$ and again we find $$ (A-I)^{-1} = \sum_{m=0}^{n-1} A^m \frac{\sum_{k=m+1}^n c_k}{\sum_{k=0}^n c_k} $$