The fact that two operators that commute have a common set of eigenfunctions can be seen as follows: assume that is an eigenfunction of with eigenvalue . Then since and commute, . Comparing start and end, must be an eigenfunction of with eigenvalue just like itself is. If there is no degeneracy of the eigenvalue, that must mean that equals or is at least proportional to it. That is the same as saying that is an eigenfunction of too. (In the special case that is zero, is still an eigenfunction of , with eigenvalue zero.)
If there is degeneracy, the eigenfunctions of are not unique and you can mess with them until they all do become eigenfunctions of too. That can be shown assuming that the problem has been approximated by a finite-dimensional one. Then and become matrices and the eigenfunction become eigenvectors. Consider each eigenvalue of in turn. There will be more than one eigenvector corresponding to a degenerate eigenvalue . Now by completeness, any eigenvector can be written as a combination of the eigenvectors of , and more particularly as where is a combination of the eigenvectors of with eigenvalue and a combination of the eigenvectors of with other eigenvalues.
The vectors and separately are still eigenvectors of if nonzero, since as noted above, converts eigenvectors of into eigenvectors with the same eigenvalue or zero. (For example, if was not , would have to make up the difference, and can only produce combinations of eigenvectors of that do not have eigenvalue .) Now replace the eigenvector by either or , whichever is independent of the other eigenvectors of . Doing this for all eigenvectors of you achieve that the replacement eigenvectors of are either combinations of the eigenvectors of with eigenvalue or of the other eigenvectors of . The set of new eigenvectors of that are combinations of the eigenvectors of with eigenvalue can now be taken as the replacement eigenvectors of with eigenvalue . They are also eigenvectors of . Repeat for all eigenvalues of .
Similar arguments can be used recursively to show that more generally, a set of operators that all commute have a common set of eigenvectors.
The operators do not really have to be Hermitian, just
diagonalizable: they must have a complete set of
The above derivation assumed that the problem was finite-dimensional, or discretized some way into a finite-dimensional one like you do in numerical solutions. The latter is open to some suspicion, because even the most accurate numerical approximation is never truly exact. Unfortunately, in the infinite-dimensional case the derivation gets much trickier. However, as the hydrogen atom and harmonic oscillator eigenfunction examples indicate, typical infinite systems in nature do satisfy the relationship anyway.