Simultaneous triangularization of matrices Let $\mathcal{F}=\{A_1,A_2,\ldots,A_r\}$ be a triangulable commuting family of $n\times n$ matrices (that is, each $A_i$ is triangulable and $A_iA_j=A_jA_i$ for every $i,j$). I know that $\mathcal{F}$ can be simultaneous triangularization, but what is the algorithm of finding the invertible matrix $P$ such that $P^{-1}A_iP$ is triangular? As a working example consider the matrices 
$$
A=
\begin{pmatrix}
-3 &amp; 2 &amp; -4 \\
-1 &amp; 0 &amp; -1\\
 2 &amp; -2 &amp; 3
\end{pmatrix}\qquad
B=
\begin{pmatrix}
3 &amp; -2 &amp; 2 \\
-1 &amp; 2 &amp; -1\\
-2 &amp; 2 &amp; -1
\end{pmatrix}\qquad
C=
\begin{pmatrix}
1 &amp; 0 &amp; 1 \\
1 &amp; 0 &amp; 1\\
0 &amp; 0 &amp; 0
\end{pmatrix}
$$
Here $AB=BA$, $AC=CA$ and $BC=CB$. In addition, the characteristic polynomials of these matrices are
$$
f_A(x)=x(x-1)(x+1)\\
f_B(x)=(x-2)(x-1)^2\\
f_C(x)=x^2(x-1)
$$
so each one of them is triangulable. Thanks!
Here, it is easy. Since $A$ has $3$ distinct eigenvalues and $B,C$ commute with $A$, we can deduce that $B,C$ are polynomials in $A$ and it suffices to triangularize $A$.

In the general case.

Step 1. Find a common eigenvector of $A,B,C$.

Step 2. Proceed by recurrence.

Moreover, you can choose $P$ as an orthogonal matrix.
More explicitly: find a common eigenvector v1 among your given matrices. Now complete a basis arbitrarily *u2,...,un. Let Q = [v,u1,...,un ] a matrix. Now Q-1AQ, Q-1BQ, and Q-1CQ should all have an (n-1) x (n-1) submatrix in the bottom right with only zeros in the entries of first column directly to the left of it. You then want to find a common eigenvector w among all of these submatrices. Let U =[u1,...,un]. Then v2 =Uw. With this update *Q = [v1,v2,u3,...,un ] *  where u3,...,un may have been rechosen in order to complete the space. You rinse and repeat with n-2, n-3 and so on. 

If you need more information check out the algorithm in the paper by Dubi: