I am trying to find out a way of simultaneously diagonalize multiple (larger) matrices, not necessarily by hand but with some help of Sagemath. Now I understand that for two matrices the method is as follows (correct me if I'm wrong):
For two matrices $A$ and $B$ (that commute), we can diagonalize $A$ and write it as $A=PDP^{-1}$. Then we compute $B_D=P^{-1}BP$ and compute the eigenvectors of $B_D$, say $v_1,...,v_n$. Then the vectors $w_1,...w_n$, where $w_i=Pv_i$ are eigenvectors for both $A$ and $B$.
But what if we now have a third matrix $C$, that commutes with both $A$ and $B$. Can we generalize the above procedure for multiple matrices?
That you wrote is false (take $B=I$).
Let $A_1,\cdots,A_k\in M_n$ be pairwise commuting and diagonalizable matrices with minimal polynomials $(\mu_{A_i})$. We want to simultaneously diagonalize these matrices; denote this problem by $\mathcal{P}_{n,k}$.
Step 1. For example, we assume that $A_1$ realizes the maximum of the degrees of the $\mu_{A_i}$. We diagonalize $A_1$ in the form $P^{-1}A_1P=diag(\lambda_1I_{n_1},\cdots,\lambda_pI_{n_p})$ where the $(\lambda_i)$ are distinct and $n_1+\cdots+n_p=n$. Then the $P^{-1}A_iP,i>1$ are in the form $diag(A_{i,1},\cdots,A_{i,p})$ where, for every $u\leq p$, $A_{i,u}\in M_{n_u}$ and, for every $i\not= j$, $A_{i,u}A_{j,u}=A_{j,u}A_{i,u}$.
We reduced $\mathcal{P}_{n,k}$ to the problems $\mathcal{P}_{n_1,k-1},\cdots,\mathcal{P}_{n_p,k-1}$.
Step 2. Proceed by recurrence.