Angular momentum commutation relations The operator $L^2$ commutes with each of the operators $L_x$, $L_y$ and $L_z$, yet $L_x$, $L_y$ and $L_z$ do not commute with each other.
From linear algebra, we know that if two hermitian operators commute, they admit complete sets of common/simultaneous eigenfunctions. The way I understand this statement is that the eigenfunctions of both operators are the same. So, if that were the case, that would mean that $L_x$ has the same eigenfunctions as $L^2$. The same goes for $L_y$ and $L_z$. That would mean that $L_x$, $L_y$ and $L_z$ all have the same eigenfunctions, which doesn't seem to be true since they do not commute with each other. How is this resolved?
 A: Your formulation is correct:

From linear algebra, we know that if two hermitian operators commute, they admit complete sets of common/simultaneous eigenfunctions.

However, if two hermitian operators commute, it's not true that every set of eigenfunctions for either of them will be one of these sets of common eigenfunctions.
Thus, $L^2$ has shared eigenbases with each of $L_x$, $L_y$ and $L_z$, but those are different eigenbases. The only state that's common to all three is the spherically-symmetric $L^2=0$ state.
A: Be very careful about what the theorem is saying! You posted in a comment:

If two Hermitian operators, A and B, commute and if A has no degenerate eigenvalue, then each eigenvector of A is also an eigenvector of B.

This statement is not symmetric in A and B! How does it apply to our situation? $L^2$ and $L_z$ commute, and $L_z$ has no degenerate eigenvalues. The statement above means that every eigenvector of $L_z$ is also an eigenfunction of $L^2$. No problem there.
But $L^2$ does have degenerate eigenvalues, and so there are some eigenvectors of $L^2$ that are not eigenvalues of $L_z$. (In particular, the eigenvectors of $L_y$ and $L_x$ are eigenvectors of $L^2$ but not of $L_z$.)
A: Commuting operators do not necessarily share ALL eigenstates, just some set.
An eigenstate shared by $L^2$ and $L_x$ will not be the same as that shared by $L^2$ and $L_y$
A: If two observables $A$ and $B$ commute, i.e. if $[A,B]=0$, then there exists a common eigenbasis. In other words, there is a basis $\{|\phi_n\rangle\}_n$ for which
$$A|\phi_n\rangle= a_n\, |\phi_n\rangle \quad\text{and}\quad  B|\phi_n\rangle= b_n\, |\phi_n\rangle \quad.$$
Now consider the case where $A$ also commutes with another observable $C$. Then this does not imply that $C|\phi_n\rangle= c_n\, |\phi_n\rangle$ for all $n$: The basis $\{|\phi_n\rangle\}_n$ is, in general, not an eigenbasis of $C$.
The fact that $L^2$ commutes with all $L_x$, $L_y$ and $L_z$ does not imply that the e.g. $L_x$ and $L_y$ commute. Indeed, as you pointed out, they do not commute and hence do not share a common eigenbasis, although each of them shares a common eigenbasis with $L^2$.
A: Since there's probably a mathematical development of this property in your textbook/notes, I'm guessing you want an intuitive approach on it. If this doesn't help you I'll develop the mathematics:
Imagine $L^2$ has having three properties, let's call them colors. $L^2$ is green, red, and blue. $L_x , L_y , L_z$ are respectively just green, red and blue. If two operators have the same color, they'll commute, and if they don't have any common color, they won't commute. You can easily see the commutation relationships you describe above are true for this "colorful" approach.
In fact, think of this funny experiment. The identity matrix surely commutes with any matrix you can think of, let's say, A and B. If it worked as you think, this would mean that any arbitrary pair of matrices A and B must commute!
A: This is possible precisely because $L^2$ is degenerate: for an eigenvalue $l(l + 1)$, it has an eigenspace of dimension $2l + 1$ (i.e., it has this many linearly independent eigenstates). The choice of basis in a degenerate eigenspace is not unique -- thus explaining how it can be that the eigenstates $L^2$ shares with $L_x$ are not the same as the ones $L^2$ shares with $L_y$.
It is possible to simultaneously diagonalize $L^2$ and $L_x$, or $L^2$ and $L_y$, or $L^2$ and $L_z$. But the most you can do for $L^2$, $L_x$, $L_y$, and $L_z$ together is simultaneous block diagonalization: $L_x$, $L_y$, and $L_z$ are nonzero only in each $(2l + 1) \times (2l + 1)$ block corresponding to each $L^2$ eigenvalue. And $L^2$ commutes with all of them because it is proportional to the identity matrix in each block.
