Invariants of co-diagonalizability in real symmetric matrices This question has been mentionned to me by U. Frisch. He wanders whether it has ever been considered by algebraists.
In the vector space ${\bf Sym}_n({\mathbb R})$, two elements commutte to each other if and only if they are diagonalisable in the same orthonormal basis. If we stay away form the singular set $\Sigma$ formed by matrices having a multiple eigenvalue (the zero set of the discriminant, it is a codimension-$2$ subset), then the co-diagonalisability is an equivalence relation, which defines equivalence classes. These classes are algebraic subsets.

What is the simplest set of rational functions $\Phi:{\bf Sym}_n({\mathbb R})\rightarrow{\mathbb R}^N$ such that these classes are level sets of $\Phi$?

An equivalence class is a sub-manifold of codimension $n(n-1)/2$. We thus expect that locally, $N=n(n-1)/2$ is enough. This might or might not be globally true.
The case $n=2$ is very easy, with $N=1$ and 
$$\Phi(S)=\frac{s_{22}-s_{11}}{s_{12}}.$$
 A: Consider the orthogonal group $G=O(n)$ acting on the space of real $n\times n$ symmetric matrices $V=Sym(n,\mathbf R)$ by conjugation, namely, $g\in O(n)$ acts on $A$ by $g\cdot A = gAg^{-1}=gAg^t$. The orbits of this action are the usual conjugacy classes of matrices.
If we equip $V$ with the trace inner product $\langle A,B\rangle =trace(AB)$
then this action becomes isometric. 
Suppose now $A$ is a regular matrix, i.e. it has pairwise different eigenvalues. Then:

The 
  co-diagonalizability class of $A$ consists precisely of the matrices lying in the normal space $N_A(G\cdot A)$ to the orbit through $A$. 

(Note that $A$ itself lies in this normal space, so the affine normal space coincides with the normal space). In fact, 
this is clear if $A$ is a diagonal matrix, since in that case the centralizer of $A$ consists of the subspace $\Sigma$ of all diagonal matrices, and $\Sigma$ is also the normal space $N_A(G\cdot A)$ as is easily checked. For a general regular matrix $A$,  there exists $g\in O(n)$ such that $gA\in\Sigma$ and then $g^{-1}\Sigma$ is both the co-diagonalizability 
class of $A$ and the normal space $N_A(G\cdot A)$.  
Suppose $A\in V$ is regular. Since $N_A(G\cdot A)$ is an equivalence class, it follows that for all regular $B\in N_A(G\cdot A)$ we also have $N_B(G\cdot B)=N_A(G\cdot A)$.
In other words:

Along the normal space $N_A(G\cdot A)$, all orbits through regular matrices have the same tangent space. 

The property above means that the action of $G$ on $V$ is an example of a polar representation. Namely, an orthogonal representation of a compact Lie group $G$ on a 
(finite-dimensional) real vector space $V$ is called polar if there exists a subspace $\Sigma$ that meets all the orbits, and meets always orthogonally. I think now it is more natural to solve the question in the context of polar representations. 
Let $(G,V)$ be a polar representation. Fix a regular point $v_0\in V$ (in this context we may assume that a regular point is just one through which the orbit has maximal dimension). Then the normal space $\Sigma=N_{v_0}(Gv_0)$ is a section. The other sections are all of the form $g\Sigma$, for some $g\in G$. Two different sections can only meet at a singular point, so that the regular points of all possible sections foliate the set of 
regular points. We want a rational map $\Phi:V\to \mathbf R^p$ whose level sets are the sections. 
For regular $v\in V$, we study the intersection of the normal space $N_v(Gv)$ with the 
affine tangent space $v_0+T_{v_0}(Gv_0)$. Note that they have complementary dimensions in $V$. We claim that for $v$ lying in the complement of a proper real algebraic set of $V$,
that intersection consists of a unique point, which we call $v_0+\Phi(v)$. This will 
define a rational map $\Phi:V\to T_{v_0}(Gv_0)\cong\mathbf R^p$, where $p$ is the dimension of maximal dimensional $G$-orbits. 
In fact, the linear system that defines $\Phi(v)$ is 
$\langle \Phi(v),X\cdot v\rangle = - \langle v_0,X\cdot v\rangle$ for $X\in\mathfrak g$.
Let $X_1,\ldots,X_m$ be a basis of the Lie algebra 
$\mathfrak g$ such that $X_{p+1},\ldots,X_m$ lies in the isotropy algebra at $v_0$. Then $\Phi(v)=\sum_{j=1}^px_jX_j\cdot v_0$ for some real numbers $x_1,\ldots,x_p$ depending on $v$, and the above system has 
the explicit form

$\sum_{j=1}^px_j\langle X_j\cdot v_0,X_i\cdot v\rangle=\langle X_i\cdot v_0,v\rangle$ for 
  $i=1,\ldots,m$.

The case $m=p$, which includes the case of symmetric matrices, is simpler in that the system has the same number of equations as variables and one can solve it uniquely whenever $v$ satisfies 
$\det(\langle X_j\cdot v_0,X_i\cdot v\rangle)\neq0$,
thus yielding a rational map $\Phi:V\to\mathbf R^p$,
$v\mapsto (x_1,\ldots,x_p)$. In the case $m>p$, one needs to analyse more carefully the system.   
The construction depends on the choices of a regular point $v_0\in V$ and a basis $X_1,\ldots,X_m$ of $\mathfrak g$, but presumably not in an essential way.  
Back to the case of $n\times n$ symmetric matrices, one can write $\Phi$ explicitly, but the formulae for general $n$ seem not to tell much. The case $n=2$ is very simple, as has been pointed out by the OP: taking 
$v_0$ to be the diagonal matrix with entries $(1,-1)$ we get $\Phi(a_{ij})=2a_{12}/(a_{11}-a_{22})$.  
