Matrices commute if and only if they share a common basis of eigenvectors? I've come across a paper that mentions the fact that matrices commute if and only if they share a common basis of eigenvectors. Where can I find a proof of this statement?

An elementary argument.
Summary: show that each eigenspace of $A$ has a basis such that each basis vector is contained in one of the eigenspace of $B$. This basis is then the simultaneous common basis we are looking for.
Suppose $A,B$ are both diagonalizable and they commute.
Now let $E_{\lambda_i}$ be eigenspaces of $A$ for each distinct eigenvalue $\lambda_i$ of $A$.
Now let $F_{s_i}$ be eigenspaces of $B$ for each distinct eigenvalue $s_i$ of $B$.
Now I claim that $E_{\lambda_i}$ (of say dimension $m$) has a basis $v_1^i,...,v_m^i\in E_{\lambda_i}$ such that each $v_r^i$ is in one of $B$'s engenspace $F_{s_j}$--this would imply these $v_r^i$ are eigenvectors of $B$ and $A$ simultaneously. Apply this to all eigenspaces $E_{\lambda_i}, i=1,...,n$. The collection of all $v_r^i$ then becomes a common basis for $A$ and $B$ as required.
To show this claim, first pick arbitrary basis $w_1,...,w_m$ of $E_{\lambda_i}$. Each $w_i$ can be written as sum of vectors where each vector is in one of $B$'s engenspace $F_{s_j}$. This is a subtle point so let me repeat: for each $i=1,...,m,$ $w_i=z_1^i+...+z_{l_i}^i, l_i\le m$ and $z_k^i\in F_{s_j}$ for some $j$. This is trivially true because direct sum of $B$'s engenspaces is the entire space.
Now we make a second claim that all $z_k^i\in E_{\lambda_i}$. Then the collection of all $z_k^i$ span $E_{\lambda_i}$ and thus the collection can be reduced to a basis $v_1,...,v_m$ where each $v_j$ is contained in $E_{\lambda_i}$ as required by the first claim.
Note that $B$ is invariant to $E_{\lambda_i}$ since $A,B$ commute. The second claim follows from: $\sum_{i=1}^N z_i \in S$ where $z_i$ are eigenvectors of distinct eigenvalues of $B$ and $S$ is a subspace to which $B$ is invariant, then $z_i\in S,\forall i$. We check this by induction on $N$. It is trivially true for $N=1$. Then suppose $Bz_1=\lambda z_1$. Since $\lambda(z_1+...+z_N)\in S$ and $B(z_1+...+z_N)\in S$, we have $B(z_1+...+z_N)-\lambda(z_1+...+z_N)=a_2z_2+...+a_Nz_N\in S$ for some constant $a_i\neq 0$--the constants are non-zero because we assumed $z_i$ all have distinct eigenvalues. Then apply inductive hypothesis $z_2,...,z_N\in S$. This would imply $z_1\in S$ as well. This finishes the proof.
