Getting a simultaneous basis of diagonalization for a list of diagonalizable matrices that
pairwise commute is just a matter of getting a diagonalization basis for any one of them!

This can be obtained with `jordan_form`

with the argument `transformation=True`

.

For example, starting from these two matrices `a`

and `b`

:

```
sage: a = matrix(QQ, [[1, 1, 1], [-2, 3, 2], [0, 1, 2]])
sage: b = matrix(QQ, [[-1, 1, 1], [1, 1, -1], [-3, 1, 3]])
```

Check that they commute:

```
sage: a * b
[-3 3 3]
[-1 3 1]
[-5 3 5]
sage: b * a
[-3 3 3]
[-1 3 1]
[-5 3 5]
sage: a * b == b * a
True
```

If we diagonalize them separately:

```
sage: da, pa = a.jordan_form(transformation=True)
sage: db, pb = b.jordan_form(transformation=True)
```

we can notice that the transformation matrices `pa`

and `pb`

have the same columns, permuted:

```
sage: pa
[ 1 1 0]
[ 1 0 1]
[ 1 1 -1]
sage: pb
[ 0 1 1]
[ 1 1 0]
[-1 1 1]
```

Obtaining the simultaneous diagonal form is just a matter of reordering the columns of:

```
sage: da
[3|0|0]
[-+-+-]
[0|2|0]
[-+-+-]
[0|0|1]
sage: db
[2|0|0]
[-+-+-]
[0|1|0]
[-+-+-]
[0|0|0]
```

For instance, use the transformation matrix `pa`

to get simultaneous diagonalizations:

```
sage: qa = ~pa
sage: qa * a * pa
[3 0 0]
[0 2 0]
[0 0 1]
sage: qa * b * pa
[1 0 0]
[0 0 0]
[0 0 2]
```

could you write the code?

Please give us at least two of the many commuting matrices that can be diagonalized (simultaneously).

Let me join the club of asksage junkies : we need your code to understand your problem and (hopefully) provide a solution.

Please provide an example of matrices for the computation you are asking about.

Ready-to-copy-paste examples make it way easier to explore a question, thereby increasing chances of an answer.