Implicit finite element method implemented on 2D heat equation diverge at small time stem size? I applied implicit FEM on the problem $$ u_t -\Delta u = f$$ using quadratic Lagrange polynomial as a basis. The system of linear equation is
$$ (M + \Delta t A)u^{n+1} = Mu^{n} + \Delta t f^n, $$ where $M_{ij} = (\phi_i,\phi_j)$ and $A_{ij} = (\nabla \phi_i,\nabla \phi_j) $. The solution and rate of convergence is correspond to theorem but the solution is blow-up in the small $\Delta t$. I think it should be unconditionally stable but the result betrayed me. I write the problem here to persuade someone who read the question to confirm the unconditionally stable condition or the result is right.
 A: This scheme is unconditionally stable in the energy norm $\|u\|_A^2 = (Au, u)$
$$
M \frac{u^{n+1} - u^n}{\Delta t} + A u^{n+1} = f^n\\
\left(M + \frac{\Delta t}{2} A\right) \frac{u^{n+1} - u^n}{\Delta t} 
- A \frac{u^{n+1} - u^n}{2} 
+ A u^{n+1} = f^n\\
\left(M + \frac{\Delta t}{2} A\right) \frac{u^{n+1} - u^n}{\Delta t} 
+ A \frac{u^{n+1} + u^n}{2} = f^n\\
$$
Now multiply both sides with $u^{n+1} - u^n$ to get
$$
\frac{1}{\Delta t} \|u^{n+1} - u^n\|_{M + \frac{\Delta t}{2}A}^2
+ \frac{1}{2} \left(\|u^{n+1}\|^2_A - \|u^n\|^2_A\right) = (f^n, u^{n+1} - u^n).
$$
$$
\|u^{n+1}\|^2_A = \|u^n\|^2_A 
- \frac{1}{\Delta t} \|u^{n+1} - u^n\|_{M + \frac{\Delta t}{2}A}^2
+ \underbrace{2 (f^n, u^{n+1} - u^n)}_{O(\Delta t)}.
$$
This shows that for $f = 0$ the $\|u^n\|_A^2$ not increases in time. And if $f$ is not zero then $\|u^n\|_A^2$ increases as much as $n \cdot O(\Delta t)$ which is limited.
The proof for $f \neq 0$ is sloppy (the $n \cdot O(\Delta t)$ part), but for the case when $f = 0$ the proof is solid and could be checked in practice directly by evaluating $\|u^n\|_A^2 = (Au^n, u^n)$ on each time step. The sequence of $\|u^n\|_A^2$ must be non increasing.
