Also known as the Backward Euler Finite Difference Method.
In the Explicit Euler exposition just prior, it is worth remembering that the second dervative with respect to $x$ is taken at time $t$. Recall that we started off with
$$ \dfrac{u(x,t + \Delta t ) \; – \; u(x,t)}{\Delta t} \simeq \dfrac{\partial{u}}{\partial{t}} = \dfrac{\partial^2{u}}{\partial{x^2}} \simeq \dfrac{u(x + \Delta x, t )\; – \; 2u(x,t) \; + \; u(x – \Delta x,t)}{(\Delta x)^2}$$
this gave us an explicit solution for the temperature at $\Delta t = 1$ or after a single time step
$$U_j^{i+1} = U_j^i + \mu (U_{j+1}^i \; – \; 2U_j^i + U_{j-1}^i) \qquad 1 \leqslant j \leqslant M-1$$
If we make a small amendment to the explicit equation by taking the second partial derivative at $t+\Delta t$
\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}
\end{align}
again where $\mu = \dfrac{k}{h^2}$.
We obtain the Implicit Euler solution and the theorem we end up with is much nicer.
Theorem
Implicit Euler is stable for ALL $\mu$.
The Matrix Form of Implicit Euler
The equation for $U_j^{i}$ may be expressed in matrix form as
\begin{align}
\mathbf{U^{i}} &= T_I\mathbf{U^{i+1}} \\ \\
T_I &=
\left(\begin{array}{ccccc} 1+2\mu & -\mu & & & \\ -\mu & 1+2\mu & -\mu & \\ & \ddots & \ddots & \ddots & \\ & & -\mu & 1+2\mu & -\mu \\ & & & & 1+2\mu\end{array}\right) \quad \in \mathbb{R}^{(M-1) \times (M -1)}\\ \\
\mathbf{U^{i+1}} &= \left(\begin{array}{c}U^{i+1}_1 \\ \vdots \\ U^{i+1}_{(M-1)}\end{array}\right) \qquad \in \mathbb{R}^{(M-1)}
\end{align}
And the matrix $\mathbf{U^i}$ is of order $M$ seconds and not order $M^3$.
Powers of $T_I$ do not blow up like powers of the explicit Euler $T$ matrix do. To prove this we need to look at the eigen values of the matrices. If $T_I$ were a general matrix, the operation above would take a multiple of $M^3$ operations in terms of time as solving linear equations is generally cubic in the number of variables. However, for a tridiagonal matrix the number of operations is linear in $M$ and therefore greatly reduced.
Next, we explore the Stability Properties of Implicit Euler