I have shown that a finite difference method satisfies
$$\underline{u}^{n+1}=((1+6\mu)\mathbb{I}-36\mu A^{-1})\underline{u}^n$$
I don't think that the rest of the question is necessary but it is a finite difference method for the diffusion equation on $t>0, 0<x<1$ with boundary conditions $u(0,t)=0=u(1,t)$ and $u(x,0)=u_0(x)$.
The last part of the question asks to show that the scheme is stable if $\mu\leq 1/6$. The $A$ matrix was a tridiagonal symmetrix matrix with the numbers $1,4,1$ down the diagonal. It is said that the eigenvalues for this matrix are
$\lambda_i=4+2\cos\left( \frac{i\pi}{M+1} \right)$ for $i=1,\cdots,M$ if it is an $M\times M$ matrix.
Attempt
If I can show that the matrix
$$\left[(1+6\mu)\mathbb{I}-36\mu A^{-1}\right]$$
Is bounded with respect to some norm then the finite difference method is stable.
So I think that the norm will be the spectral norm since we have been given eigenvalues
If the eigenvalues are as above then take the reciprocal. So the eigenvalues of $A^{-1}$ is $\frac{1}{4+2\cos(\cdots)}$ which has a maximum value of 1/2.
So now using this I will try to show that
$$||B||=||(1+6\mu)\mathbb{I}-36\mu A^{-1}||=\sqrt{\rho(B^TB)}\leq C$$
I get stuck in evaluating this. My guess would be that the matrix is equal to its transpose so it is just:
$$||B||=\sqrt{\rho({B^2})} $$
My linear algebra is not very good I don't know if you can do this last part. Maybe take
$$(1+6\mu)^2\mathbb{I}-(1+6\mu)A^{-1}+36^2\mu^2A^{-1}A^{-1}$$
and if it were to act on some eigenvector then the eigenvalue would be
$$(1+6\mu)^2-\frac{1}{2}(1+6\mu)+\frac{36^2 \mu^2}{4}=\frac{1}{2}+9 \mu +360 \mu ^2$$
This doesn't seem to shed any light on it.
On the other hand, if I just take
$$|\rho(B)|=|1+6\mu-36\mu\frac{1}{2}|=|1-12\mu|$$
Then I for $0<\mu \leq 1/6$
I have $||B|| \leq 1$?
You're right that the eigenvalues of $A^{-1}$ are $\lambda^{-1}$ for each eigenvalue $\lambda$ of $A$. (The notation is suggestive for a reason.) Furthermore, the identity matrix commutes with every matrix; so the eigenvalues of
$qI-rA^{-1}$
are just $q-r/\lambda$ for each eigenvalue $\lambda$ of $A$. Using this, you can compute the eigenvalues of $(1+6\mu)I-36\mu A^{-1}$.
You also know that the matrix $A$ is symmetric, so its inverse is also symmetric. For a symmetric matrix, the 2-norm $\|B\|_2 = \sqrt{\rho(B^\top B)}$ is just the absolute value of the eigenvalue of largest magnitude:
$\|B\|_2 = \max\{|\lambda| : \lambda \textrm{ is an eigenvalue of }B\}$.
I think you've more or less done this at the end of your answer without knowing the justification for why you can do it. You may also want to look at the smallest eigenvalue of $A^{-1}$.
Edit: There are actually plenty of definitions of stability, but I think the one you're looking for is that the solution's norm is bounded as $n\to\infty$. For that, you need that $\|B\| \le 1$ in some norm, where $B$ is the matrix corresponding to the iteration. Provided that you have a consistent discretization of your PDE, the Lax-Richtmeyer Theorem guarantees that, if the discretization is also stable, it will converge to the true solution as the mesh spacing $h$ goes to 0.