Substituting matrices in polynomials

This issue is inspired by [this](https://math.stackexchange.com/questions/4250154/cayley-hamilton-in-macaulay2) Stack Exchange question that I randomly ran into.

Would be great if the following worked:
```m2
S = QQ[t]
m = random(ZZ^2, ZZ^2)
f = det(m - t*m^0) -- characteristic polynomial
-- By Cayley-Hamilton theorem we should have:
0 == f[m] -- error: expected a matrix with 1 row
0 == sub(f, t => m) -- error: expected substitution values to be ring elements or numbers
```
The easiest workaround I can think of is still too complicated:
```m2
0 == sum(pairs reverse flatten entries last coefficients(f, Variables => t), (i, c) -> m^i * c)
```
But more generally, if I have a polynomial in some number of variables, it shouldn't be too difficult to substitute matrices (or even ring maps) for the variables or return an error if the matrices are not compatible in some way.

Related: #2863, #2792

Start by reproducing the two failing examples in the issue and read related issues #2863 and #2792. Define the supported substitution behavior and compatibility errors for matrix-valued variables, then verify that the Cayley-Hamilton example and the broader multivariable case work as requested.

Written by the indexing model from the issue text.