Introduction over unbounded domain
Consider the forward time shift $\mathsf{z}$ acting on a discrete function of time (sequence) $f=(f^n)_{n\in\mathbb{N}}$ as $(\mathsf{z} f)^{n} = f^{n+1}$. Also consider the backward space shift $\mathsf{s}$ acting on a discrete function of space by $f=(f_{j})_{j\in\mathbb{Z}}$ as $(\mathsf{s} f)_{j} = f_{j-1}$.
We consider the following vectorial iteration $$ \mathbf{y}^{n+1} = \mathbf{A} \mathbf{y}^{n}, \qquad \text{where} \qquad \mathbf{A} \in \mathcal{M}_{q}(\mathbb{R}[\mathsf{s}, \mathsf{s}^{-1}]), $$ where $\mathbb{R}[\mathsf{s}, \mathsf{s}^{-1}]$ is the ring of Laurent polynomials in the indeterminate $\mathsf{s}$ over real numbers.
To make things concrete, we will use $$ \mathbf{A} = \begin{bmatrix} \tfrac{1}{2}(\mathsf{s} + \mathsf{s}^{-1}) & -\tfrac{1}{2}(\mathsf{s} - \mathsf{s}^{-1}) \\ \tfrac{1}{2}(\mathsf{s} - \mathsf{s}^{-1}) & -\tfrac{1}{2}(\mathsf{s} + \mathsf{s}^{-1}) \end{bmatrix}. $$
If the aim is to rewrite, far from the initial time $n = 0$, the dynamics of $y_1$ (the first component of $\mathbf{y}$) only as a function of itself, one has to rewrite the system as $(\mathsf{z}\mathbf{I}-\mathbf{A})\mathbf{y}^n = \mathbf{0}$, multiply it by the adjugate $\text{adj}(\mathsf{z}\mathbf{I}-\mathbf{A})$ and invoking the well-known formula for the adjoint $\text{adj}(\mathbf{C})\mathbf{C} = \text{det}(\mathbf{C}) \mathbf{I}$, which is indeed equivalent to the Cayley-Hamilton theorem. This gives $\text{det}(\mathsf{z}\mathbf{I}-\mathbf{A})\mathbf{y}^n = \mathbf{0}$. Selecting the first row yields the desired result: $$ \text{det}(\mathsf{z}\mathbf{I}-\mathbf{A})y_1^n = 0. $$
On the previous example $$ \text{det}(\mathsf{z}\mathbf{I}-\mathbf{A}) = \mathsf{z}^2 - \text{tr}(\mathbf{A}) \mathsf{z} + \text{det}(\mathbf{A}) = \mathsf{z}^2 + 1. $$ This corresponds to the three-level recursive iteration $$ y_{1, j}^{n+2} + y_{1, j}^{n} = 0. \qquad \qquad (A) $$
This process using the characteristic polynomial of $\mathbf{A}$ provides an "automatic" way of doing the substitutions one can do by hand.
Write the equations of the problem: \begin{align} y_{1,j}^{n+1} = \tfrac{1}{2}(y_{1,j-1}^{n}+y_{1,j+1}^{n}) - \tfrac{1}{2}(y_{2,j-1}^{n}-y_{2,j+1}^{n}), \qquad \qquad (1)\\ y_{2,j}^{n+1} = \tfrac{1}{2}(y_{1,j-1}^{n}-y_{1,j+1}^{n}) - \tfrac{1}{2}(y_{2,j-1}^{n}+y_{2,j+1}^{n}). \qquad \qquad (2) \end{align} Consider (2) at the previous time step shifting back and forth in space, in practice $(2)^{n-1}_{j-1}-(2)^{n-1}_{j+1}$ $$ y_{2,j-1}^{n} - y_{2,j+1}^{n}= \tfrac{1}{2}(y_{1,j-2}^{n-1}-2y_{1,j}^{n-1}+y_{1,j+2}^{n-1}) - \tfrac{1}{2}(y_{2,j-2}^{n-1}-y_{2,j+2}^{n-1}). \qquad \qquad (3) $$ We still have to get rid of $y_2$ in the previous equation. To this end, we take $(1)^{n-1}_{j-1}+(1)^{n-1}_{j+1}$:
$$ y_{1,j-1}^{n} + y_{1,j+1}^{n} = \tfrac{1}{2}(y_{1,j-2}^{n-1}+2 y_{1,j}^{n-1}+y_{1,j+2}^{n-1}) - \tfrac{1}{2}(y_{2,j-2}^{n-1}-y_{2,j+2}^{n-1}). $$ We isolate the last term $$ - \tfrac{1}{2}(y_{2,j-2}^{n-1}-y_{2,j+2}^{n-1}) = y_{1,j-1}^{n} + y_{1,j+1}^{n} - \tfrac{1}{2}(y_{1,j-2}^{n-1}+2 y_{1,j}^{n-1}+y_{1,j+2}^{n-1}). \qquad \qquad (4) $$ We plug (4) into (3): $$ y_{2,j-1}^{n} - y_{2,j+1}^{n}= -2y_{1,j}^{n-1} + y_{1,j-1}^{n} + y_{1,j+1}^{n}. \qquad \qquad (5) $$ Inserting (5) into (1), we obtain: $$ y_{1,j}^{n+1} = y_{1,j}^{n-1} , $$ which is exactly (A) shifted at the previous time step, which is harmless since the system is time-invariant.
The issue
Now we introduce a space boundary at $j=0$ and consider a semi-infinite domain on the right. The standard iteration (1) and (2) will be used for $j \geq 1$, whereas we introduce a different rule for $j=0$.
A first choice we consider is: \begin{align} y_{1, 0}^{n+1} = \tfrac{1}{2}(y_{1, 0}^n+y_{1, 1}^n) - \tfrac{1}{2}(y_{2, 0}^n-y_{2, 1}^n), \qquad \qquad (6) \\ y_{2, 0}^{n+1} = \tfrac{1}{2}(y_{1, 0}^n-y_{1, 1}^n) - \tfrac{1}{2}(y_{2, 0}^n+y_{2, 1}^n). \qquad \qquad (7) \end{align} Again, we want to eliminate $y_2$ from (6). Taking $(7)^{n-1}-(2)^{n-1}_{j=1}$ gives $$ y_{2, 0}^{n} - y_{2,1}^{n} = -\tfrac{1}{2}(y_{1, 1}^{n-1}-y_{1,2}^{n-1}) - \tfrac{1}{2}(y_{2, 1}^{n-1}-y_{2,2}^{n-1}). \qquad \qquad (8) $$ The last term in the previous equation still need to be eliminated. Considering $(1)^{n-1}_{j=1} - (6)^{n-1}$
$$ y_{1,1}^{n} - y_{1, 0}^{n} = -\tfrac{1}{2}(y_{1, 1}^{n-1} - y_{1, 2}^{n-1}) - \tfrac{1}{2}(y_{2, 1}^{n-1}-y_{2,2}^{n-1}). $$ Isolating the desired term $$ - \tfrac{1}{2}(y_{2, 1}^{n-1}-y_{2,2}^{n-1}) = y_{1,1}^{n} - y_{1, 0}^{n} +\tfrac{1}{2}(y_{1, 1}^{n-1} - y_{1, 2}^{n-1}). \qquad \qquad (9) $$ Inserting (9) in (8) gives: $$ y_{2, 0}^{n} - y_{2,1}^{n} = - (y_{1, 0}^{n} - y_{1,1}^{n}). \qquad \qquad (10) $$ With (10) back into (6): $$ y_{1, 0}^{n+1} = y_{1, 0}^n. $$ Observe that this has been possible by the joint action of the scheme inside the domain ($j\geq 1$) and the boundary scheme ($j=0$).
First question: how is it possible to do this procedure to eliminate all the variables except the first one at the boundary automatically? In other words, there is a result analogous to the Cayley-Hamilton theorem in this case? My thought: since now the system is no longer space invariant (the matrix A depends on the point), it cannot be written using a commutative ring. We especially loose the commutative character.
A second choice we consider is: \begin{align} y_{1, 0}^{n+1} = \tfrac{1}{2}(2 y_{1, 0}^n+y_{1, 1}^n) - \tfrac{1}{2}(2 y_{2, 0}^n-y_{2, 1}^n), \qquad \qquad (11) \\ y_{2, 0}^{n+1} = \tfrac{1}{2}(2 y_{1, 0}^n-y_{1, 1}^n) - \tfrac{1}{2}(2 y_{2, 0}^n+y_{2, 1}^n). \qquad \qquad (12) \end{align} Taking $2 \times (12)^{n-1} - (2)^{n-1}_{j=1}$ $$ 2y_{2, 0}^{n} - y_{2,1}^{n} = \tfrac{1}{2} (3y_{1, 0}^{n-1}-2y_{1, 1}^{n-1} + y_{1,2}^{n-1}) - \tfrac{1}{2}(3 y_{2, 0}^{n-1}+2 y_{2, 1}^{n-1} - y_{2,2}^{n-1}) $$ Again, we try to get rid of the last term. Taking $\delta \times (11)^{n-1} + \gamma \times (1)^{n-1}_{j=1}$ provides $$ \text{something} = \text{something} - \tfrac{1}{2} ((2\delta + \gamma)y_{2, 0}^{n-1} - \delta y_{2, 1}^{n-1} - \gamma y_{2, 2}^{n-1}). $$ Hence we want to take $\delta = -2$ and $\gamma = 1$. Therefore $2\delta + \gamma = -3 \neq 3$, which indicates that it is not possible to get rid of $y_2$!
Remark: it seems that there are cases where the elimination of all the variables except the first one at the boundary is not possible. This seems to suggest that there is no analoguous result to the Cayley-Hamilton theorem. Second question: how to recognize the cases where the elimination procedure is not possible?