I know that if $A$ and $B$ are two Hermitian matrices, then $A B= B A$ if and only if their eigenspaces coincide [1]. In order to apply this test one need to compute eigenvectors of both $A$ and $B$ and then determine if they are linearly independent, which is computationally expensive.
The other potentially applicable test for commuting is simultaneous triangularisability. Two matrices are simultaneously triangularisable if there exist basis in which they both take upper triangular form. I do not know how to use this property in order to numerically establish commuting property of two matrices. However, I believe that direct implementation of such a test will be very computationally expensive.
Are there any other criteria I could encode in MATLAB in order to determine whether two (Hermitian, real symmetric) matrices commute?