Proposition. Let $R$ be a commutative ring, $n\in\mathbb N$ a natural number and let $A\in\operatorname{Mat}_{n\times n}(R)$ be an $n\times n$-matrix with coefficients in $R$. Let $P_A\in R[t]$ be the characteristic polynomial of the square matrix $A$. Suppose $f\in R[t]$ is a polynomial such that $f(A)=0$, then there exists a nature number $N\in\mathbb N$ such that $P_A$ divides $f^N$ in $R[t]$.
I wonder whether we can prove this without axiom of choice. Let me sketch a proof using axiom of choice:
We start with the special case when $R$ is a field. Consider the splitting field $K$ of the characteristic polynomial. The matrix $A$ is similar to a upper triangular matrix $B\in\operatorname{Mat}_{n\times n}(K)$. Then we have $f(B)=0$, which implies that for each eigenvalue $\lambda\in K$ of $A$, we have $f(\lambda)=0$. Consequently, $P_A$ divides $f^n$ in $K[t]$, and therefore also in $R[t]$. Note that in this special case, we don't need the axiom of choice.
Next, if $R$ is an integral domain, by passing to its fractional field $k$, we deduce that $P_A$ divides $f^N$ in $k[t]$. Since $P_A$ is monic, $P_A$ divides $f^n$ in $R[t]$.
In general, it follows from Cayley-Hamilton theorem that there exists a map of $R$-algebras $S:=R[t]/(P_A)\to\operatorname{Mat}_{n\times n}(R)$ which maps $t$ to $A$. For any prime ideal $\mathfrak p\in\operatorname{Spec}(R)$, by considering the base change $R\to R/\mathfrak p$, we deduce that $f^n=0$ in $S/\mathfrak pS$, that is, $f^n\in\mathfrak pS$ for all prime ideals $\mathfrak p\in\operatorname{Spec}(R)$. In particular, for all prime ideals $\mathfrak q\in\operatorname{Spec}(S)$, we have $f^n\in(\mathfrak q\cap R)S\subseteq\mathfrak q$. Now we deduce that $f^n$ is nilpotent (here we use the axiom of choice) in $S$, which implies the result.
My impression is that we can avoid the axiom of choice by not passing to the prime ideals. I am looking forward to such a solution.
An additional question: can we take $N=n$ in general, that is to say, $P_A$ divides $f^n$ in $R[t]$ for general commutative rings $R$?
Remark
As one user indicated in the comment, we can pass to the universal case, that is, $R=\mathbb Z[(X_{i,j})_{1\le i,j\le n},(Y_k)_{0\le k\le d}]/I$, $A=(X_{i,j})_{1\le i,j\le n}$ and $f(t)=\sum_{k=0}^dY_kt^k\in R[t]$, where $I$ is the ideal generated by the defining equations associated to $f(A)=0$ (there are $n^2$ equations, therefore this ideal is generated by $n^2$ elements).
I don't know whether $R$ is integral, or whether it is reduced. However, by applying the preceding argument, we can extract a natural number $N\in\mathbb N$ and a concrete polynomial $g\in R[t]$ such that $f^N=P_Ag$. Note that $R$ is of finite type and a theoretic elimination of axiom of choice is probable.
Finally, I find a solution. In fact, we can always pick $N=n$.
The trick is to generalize the relation that $P_A$ divides $f^n$ to an equation which is true without the restraint that $f(A)=0$.
Instead of presenting a formal argument, let me describe how we can discover it. We drop the condition that $f(A)=0$ and try to establish an equation which contains $P_A$, $f^n$ and the characteristic polynomial $P_{f(A)}$ of the matrix $f(A)$.
To do so, let's start with the special case that $R$ is a field in which the characteristic polynomial $P_A$ of $A$ splits, that is to say, there exists a tuple $\lambda_1,\dots,\lambda_n\in R$ such that $P_A(t)=(t-\lambda_1)\cdots(t-\lambda_n)$.
As in the original post, we can find an upper triangular matrix with diagonal elements being $\lambda_1,\dots,\lambda_n$ which is similar to the matrix $A$. Thanks to such a matrix, we can deduce that the characteristic polynomial $P_{f(A)}(t)$ of the matrix $f(A)$ factors as $(t-f(\lambda_1))\cdots(t-f(\lambda_n))$.
Then we have $P_{f(A)}(f(t))=(f(t)-f(\lambda_1))\cdots(f(t)-f(\lambda_n))=(t-\lambda_1)\cdots(t-\lambda_n)g(t,\lambda_1,\dots,\lambda_n)=P_A(t)g(t,\lambda_1,\dots,\lambda_n)$ where $g\in R[t,Y_1,\dots,Y_n]$ is a polynomial (of which the coefficients are in fact computed by $\mathbb Z$-polynomials in terms of coefficients of $f$). It is important to remark that $g$ is symmetric in variables $Y_1,\dots,Y_n$, therefore there exists a polynomial $h\in R[t,Z_1,\dots,Z_n]$ such that $h(t,\sigma_1(Y),\dots,\sigma_n(Y))=g(t,Y_1,\dots,Y_n)$ where $\sigma_i$ are elementary symmetric polynomials. Note that $\sigma_i(\lambda)$ are coefficients of $P_A$.
In short, we can establish an identity $P_{f(A)}(f(t))=P_A(t)H(t)$ in $R[t]$ where coefficients of $H$ are computed by $\mathbb Z$-polynomials of coefficients of $f$ and $P_A$.
Then we pass to the universal case that $R=\mathbb Z[(X_{i,j})_{1\le i,j\le n},(T_k)_{0\le k\le d}]$, $A=(X_{i,j})$ and $f(t)=\sum_k T_kt^k$. The previous arguments allows us to establish a polynomial equation $P_{f(A)}(f(t))=P_A(t)H(t)$ which holds in the algebraic closure of the field $\mathbb Q((X_{i,j}),(T_k))$ and therefore also holds in the ring $R$. By specialization we claim that this identity is also valid in all commutative rings $R$ with a matrix $A\in\operatorname{Mat}_{n\times n}(R)$ and a polynomial $f\in R[t]$.
(Somebody might challenge that the existence of the algebraic closure depends on the axiom of choice, but here in fact we don't need that: the characteristic polynomial $P_A$ is separable in the field $\mathbb Q((X_{i,j}),(T_k))$, therefore it suffices to take the splitting field of $P_A$)
In particular, if $f(A)=0$, then $P_{f(A)}(t)=t^n$ and therefore $P_{f(A)}(f(t))=f(t)^n$, the identity leads to the conclusion that $P_A$ divides $f^n$. Q.E.D.