Implicit Euler (Backward Euler):
\begin{align}U_j^{i+1} &= U_j^i + \mu (U_{j+1}^{i+1} \;- \; 2U_j^{i+1} + U_{j-1}^{i+1}) \qquad 1 \leqslant j \leqslant M-1\\ \\
\Longrightarrow \quad U_j^i &= – \mu U_{j+1}^{i+1} + (1+2\mu)U_j^{i+1} – \mu U_{j-1}^{i+1} \tag{1}
\end{align}
The characteristics of tridiagonal symmetric Toeplitz (TST) matrices hold for all $T(a,b)$ and we have an expression for the corresponding eigenvalues.
\begin{align}
T (a,b)&=
\left(\begin{array}{ccccc} a & b & & \\ b & a & b \\ & \ddots & \ddots & \ddots & \\ & & b & a & b \\ & & & & a\end{array}\right)
\quad \in \mathbb{R}^{(M-1) \times (M -1)}\\ \\
\end{align}
$$\lambda_k = a + 2b \cos \frac{\pi k}{n + 1}\; , \qquad k = 1,2,\cdots ,n $$
and this holds for $T(1+2\mu, \mu)$ may be expressed in matrix form as
\begin{align}
\mathbf{U^{i}} &= T_I(1+2\mu, \mu)\mathbf{U^{i+1}} \tag{2}
\end{align}
In general, as we saw in Stability Problems of Explicit Euler (forward Euler) any eigenvalues $\lambda_k \gt 1$ in modulus i.e outside the range $[-1,1]$ would cause exponential instability. Remember that $U_j^{i+1}$ is a numerical approximation to the PDE (such as $\dfrac{\partial{u}}{\partial{t}}$) that tells us the state at time $i+1$ of the time-space variable $\dfrac{\partial^2{u}}{\partial{x}^2}$ (or other second partial derivative) at location $jh, (h \equiv \Delta x)$ from the specified boundary. In short, we are interested in the evolution of $u(x,t)$ through time and space in a stable manner. Stability in this sense simply means the variables should not disobey the laws of physics.
Evolution of $\mathbf{U}$, the approximation of $\partial{u(x, t)}/\partial{t}$
We don’t need facts about the future to tell us about what we know now. Observe carefully,this is what the equations $(1)$ and $(2)$ seem to suggest we are doing. Meanwhile, we can and still need to determine the eigenvalues $\lambda_k$ of $T_I(1+2\mu, -\mu)\equiv T (a,b)$. So since
$$\lambda_k = a + 2b \cos \frac{\pi k}{n + 1}\; , \qquad k = 1,2,\cdots ,n $$
for $T (a,b)$, it follows that for$T_I(1+2\mu, -\mu)$,
$$\lambda_k = 1+2\mu – 2\mu \cos \frac{\pi k}{M}\; , \qquad k = 1,2,\cdots ,n $$
Clearly the function will reach its smallest value when $\cos\dfrac{\pi k}{M}=1$ such that
$$\lambda_k = 1 + 2\mu – 2\mu = 1$$
so
$$\lambda_k = 1+2\mu – 2\mu \cos \frac{\pi k}{M} \geqslant 1$$
Given what we know about explicit Euler this initially seems worrying.
What must the eigenvalues be? It all depends on the evolutionary representation of the finite difference method in question
Recall that with explicit Euler,
\begin{align}U_j^{i+1} &= \mu U_{j+1}^i + (1-2\mu)U_j^i + \mu U_{j-1}^i \qquad 1 \leqslant j \leqslant M-1\\ \\
\mathbf{U^{i+1}} &= T(1-2\mu, \mu) \mathbf{U^i}
\end{align}
we had no requirement to modify this explictly evolutionary equation as it explicitly gives future $i+1$ state based on state $i$. And the eigenvalues of $T(1-2\mu, \mu)$ between $-1$ and $1$ ensure stability.
For implicit Euler above, we need to invert $T_I(1+2\mu, -\mu)$ to obtain a forward evolution form of the equation because as it stands it is looking backward to $i$ from $i+1$ (using a future state to determine a current state). Hence $(2)$ becomes,
\begin{align}
\mathbf{U^{i+1}} = T_I(1+2\mu, \mu)^{-1} \mathbf{U^{i}} \tag{3}
\end{align}
In this form, what is required are the eigenvalues $\lambda_{k}(T^{-1})$ of $T^{-1}$. The eigenvalues of the inverse of $T_I$ are the inverse of the eigenvalues of $T_I$. That is;
$$\lambda_{k}(T^{-1}) = \dfrac{1}{\lambda_k} = \dfrac{1}{1+2\mu – 2\mu \cos \dfrac{\pi k}{M}} < 1 $$ Because the eigenvalues of $T_I(1+2\mu, \mu)$ are always $>1$, it follows that the eigenvalues $1/\lambda_k(T)$ of $T_I(1+2\mu, \mu)^{-1}$ will always be less than $1$. With implicit Euler, it is $T_I(1+2\mu, \mu)^{-1}$ that gives future $i+1$ state based on state $i$.