Solving a system of $2n-2$ equations, what is the relationship of solutions between successive n? I am trying to solve a system of $2n-2$ equations.  The first two rows of the matrix representation are always $$\begin{pmatrix}-8&1&1&0&\ldots& 0\end{pmatrix} \text{and}\begin{pmatrix}-1&-7&2&1&0&\ldots& 0\end{pmatrix}.$$
The bottom rows are always $$\begin{pmatrix}0&\ldots& 0&1&2&-6&2\end{pmatrix} \text{and}\begin{pmatrix}0&\ldots& 0&1&2&-6\end{pmatrix}.$$  The rows inbetween are $\begin{pmatrix}1&2&-6&2&1&0&\ldots& 0\end{pmatrix}$ and shifting over one to the right as you go down.  In my particular problem $x_1$ is always $1$ so that's why there's only $2n-2$ variables.
For example, for $n=4$, the matrix equation is:
$$
M_4\bf{x}=\begin{pmatrix}-8 & 1 & 1 & 0 & 0 & 0   \\
  1 &-7&  2&  1 & 0 & 0   \\
  1 & 2 &-6 & 2 & 1 & 0   \\
  0 & 1 & 2 &-6 & 2 & 1   \\
  0 & 0 & 1 & 2 &-6 & 2  \\
  0 & 0 & 0 & 1 & 2 &-6 \\
\end{pmatrix}_{6\times6}\begin{pmatrix}x_2\\x_3\\x_4\\x_5\\x_6\\x_7\end{pmatrix}=\begin{pmatrix}a\\b\\c\\d\\e\\f\end{pmatrix}
$$
For $n=5$
$$
M_5\bf{x}=\begin{pmatrix}-8 & 1 & 1 & 0 & 0 & 0 & 0 & 0 \\
  1 &-7&  2&  1 & 0 & 0 & 0 & 0 \\
  1 & 2 &-6 & 2 & 1 & 0 & 0 & 0 \\
  0 & 1 & 2 &-6 & 2 & 1 & 0 & 0 \\
  0 & 0 & 1 & 2 &-6 & 2 & 1 & 0\\
  0 & 0 & 0 & 1 & 2 &-6 & 2 & 1\\
  0 & 0 & 0 & 0 & 1 & 2 &-6 & 2\\
  0& 0 & 0 & 0 & 0 & 1 & 2 &-6\\
\end{pmatrix}_{8\times 8}\begin{pmatrix}x_2\\x_3\\x_4\\x_5\\x_6\\x_7\\x_8\\x_9\end{pmatrix}=\begin{pmatrix}a\\b\\c\\d\\e\\f\\g\\h\end{pmatrix}
$$
My particular problem is that I want to solve thousands (possibly millions) of such systems but it is computationally demanding.  You'll notice that you can step down - the $6\times6$ is simply the $8\times8$ matrix with the last two rows and columns removed.  Similarly, you can step up by extending the previous matrix with two rows and columns.
The core of my question is this: is there some relationship between the inverse of $M_n$ and $M_{n+1}$ that I can utilize to make solving more efficient?
Specifically, I am interested in $x_n$ for each $n$ if it is more efficient to only calculate that.
 A: Consider the following symmetric "pentadiagonal" Toepliz matrix :
$$T=\begin{pmatrix}-6 & 2 & 1 & 0 & 0 & 0 & 0 & 0 \\
  2 &-6&  2&  1 & 0 & 0 & 0 & 0 \\
  1 & 2 &-6 & 2 & 1 & 0 & 0 & 0 \\
  0 & \ddots&\ddots & \ddots &\ddots & \ddots & 0  & 0 \\
  0 & 0 &\ddots & \ddots &\ddots & \ddots & \ddots & 0\\
  0 & 0 & 0 & 1 & 2 &-6 & 2 & 1\\
  0 & 0 & 0 & 0 & 1 & 2 &-6 & 2\\
  0& 0 & 0 & 0 & 0 & 1 & 2 &-6\\
\end{pmatrix}_{n \times n}$$
The matrix you describe in the general case is a so-called "rank-2 perturbation" of matrix $T$, which is
$$M=T-UU^T$$
where
$$U^T=\begin{pmatrix}1&1&0&0&\cdots 0 \\ 1&0&0&0&\cdots 0 
\end{pmatrix}$$
Such a perturbation of a matrix gives rise to a similar formula giving the induced perturbation on their inverses in terms of Woodbury Matrix Identity whose general form is
$$\displaystyle{M^{-1}=\left(T+UCV\right)^{-1}=T^{-1}-T^{-1}U\left(C^{-1}+VT^{-1}U\right)^{-1}VT^{-1}}$$
which boils down in our case, taking $C=-I_n$ (identity matrix of order $n$) and $V=U^T$ :
$$M^{-1}=\left(T-UU^T\right)^{-1}=T^{-1}-\underbrace{T^{-1}U\left(-I+U^TT^{-1}U\right)^{-1}U^TT^{-1}}_{\text{perturbation term } P}$$
Dealing with inversion of matrix $T$ itself, see for example this article "An explicit formula for the inverse of a pentadiagonal
Toeplitz matrix".
Edit : In fact, even for moderate values of $n$, $T^{-1}$ is close to $M^{-1}$. Take a look at the following picture displaying $T^{-1}$ on the left and $M^{-1}$ on the right ($50 \times 50$ matrices) under a colored form with a conversion bar on the left. We see that the differences are very tiny, which is confirmed by precisely comparing the different entries (no entry in the perturbation matrix $P$ is larger than $10^{-4}$).

