Since they are all simultaneously diagonalizable, there exists a basis in which they are all diagonal. We fix this as our "standard" basis. Suppose in this basis $M_i$ has $d_{j; i}$ on as the $j$th diagonal entry. Then a linear combination $\sum c_i M_i$ is also diagonal and has $\sum c_i d_{j;i}$ as the $j$th diagonal entry.
Now what we want will follow from the two lemmas:
Lemma 1: A diagonal matrix stays diagonal in precisely those bases for which the change of basis matrix is block diagonal, with blocks corresponding to the same values on the diagonal.
Lemma 2: For "generic" $(c_1, \ldots, c_n)$ and any fixed pair $j_1, j_2$ we have $\sum c_i d_{j_1;i}=\sum c_i d_{j_2;i}$ if and only if for all individual $i$ we have $d_{j_1;i}=d_{j_2;i}$.
Indeed, the lemmas (which we prove later, making the notion of "generic" clear in the process) imply that for a generic $(c_1, \ldots, c_n)$ the only bases in which the sum is diagonal are those differing from the standard one by block-diagonal change-of-basis matrices which also keep all individual $M_i$ diagonal as well. This is exactly saying that diagonalizing the sum will diagonalize all of the $M_i$s.
Proof of lemma 1: A matrix is diagonal precisely in those bases that consist of eigenvalues. Any two such bases differ by changing bases for each eigenspace separately, meaning the change of bases matrix is block diagonal, with blocks corresponding to distinct eigenvalues. QED.
Proof of lemma 2: One way implication is clear: if individual entries are equal then the linear combinations are equal, for all $\vec{c}=(c_1, \ldots, c_n)$.
Suppose now for a given pair $j_1$ $j_2$ there exists $k$ such that $d_{j_1,k}\neq d_{j_2,k}$. Then the equation $\sum c_i d_{j_1;i}=\sum c_i d_{j_2;i}$ is of the form $c_k(d_{j_1;k}-d_{j_2;k})=\sum_{i\neq k} c_i(d_{j_1;i}-d_{j_2;i})$ and is a non-trivial linear equation. This means it is satisfied by a codimension one set of $\vec{c}$s.
Thus for all $\vec{c}$ outside of a finite collection of hyperplanes an equation $\sum c_i d_{j_1;i}=\sum c_i d_{j_2;i}$ implies $d_{j_1;i}=d_{j_2;i}$ for all $i$. Those are our "generic" $\vec{c}$s, and for them the lemma holds. QED.
Let's assume that you can show that only a diagonal matrix commutes with all diagonal matrices (this is a triviality, by comparing matrix coefficients). Now notice that matrix that commutes with all matrices must be diagonal with respect to all bases (conjugate and repeat the previous argument). This means that every vector is an eigenvector, which, in turn, obviously means that the matrix is a multiple of the identity.
Best Answer
The answer is yes, a collection of commuting diagonalisable matrices admit an basis which is an eigenbasis of each matrix. To think about why, it's better to think about linear operators and eigenspaces.
Let's say we have $A, B, C \colon V \to V$ three linear operators on a finite-dimensional vector space $V$, which pairwise commute and each is diagonalisable. For each $\lambda$, let $V(\lambda) = \{v \in V \mid Av = \lambda v\}$ be the $\lambda$-eigenspace of $A$. If $\lambda \neq \mu$ it is easy to prove that $V(\lambda) \cap V(\mu) = \{0\}$, and since $A$ is diagonalisable we have $$ V = \bigoplus_{\lambda \in K} V(\lambda),$$ where $K$ is the field you are working over (for example $\mathbb{R}$ or $\mathbb{C}$). Note that the sum above looks infinite, but there are only finitely many $\lambda$ such that $V(\lambda) \neq 0$, so it is in fact a finite sum. You should think of the operator $A$ as cutting $V$ up into pieces, each piece labelled by its eigenvalue.
The inductive step is this: since $B$ commutes with $A$, we have $B(V(\lambda)) \subseteq V(\lambda)$ since if $v \in V(\lambda)$ then $$ A(Bv) = B(Av) = B(\lambda v) = \lambda (Bv). $$ Therefore $B$ restricts to a linear operator on each $V(\lambda)$, and the same goes for $C$. So for each eigenvalue $\lambda$, we have the operators $B|_{V(\lambda)} \colon V(\lambda) \to V(\lambda)$ and $C|_{V(\lambda)} \colon V(\lambda) \to V(\lambda)$, which are commuting diagonalisable operators on the subspace $V(\lambda)$. Now applying the same logic as above, the operator $B|_{V(\lambda)}$ cuts the space $V(\lambda)$ up into pieces, each labelled by an eigenvalue of $B$. Let's name these: $$ \begin{aligned} V(\lambda, \mu) &= \{ v \in V(\lambda) \mid Bv = \mu v \} \\ &= \{v \in V \mid Av = \lambda v \text{ and } Bv = \mu v \}. \end{aligned}$$ Since $B$ is diagonalisable this sum is complete, so we have $$ V(\lambda) = \bigoplus_{\mu \in K} V(\lambda, \mu), $$ and by putting all of the $V(\lambda)$ back together we get $$ V = \bigoplus_{\lambda, \mu \in K} V(\lambda, \mu). $$ Now the fact that $C$ commutes with both $A$ and $B$ means that $C$ preserves each simultaneous eigenspace $V(\lambda, \mu)$, and so we do the same thing. It should be clear that as long as you have finitely many linear operators, you can carry out this process to completion. For our $A, B, C$ here we will get the decomposition $$V = \bigoplus_{\lambda, \mu, \nu \in K} V(\lambda, \mu,\nu),$$ where $V(\lambda, \mu, \nu)$ consists of all vectors $v$ for which $Av = \lambda v$, $B v = \mu v$, and $Cv = \nu v$. As we let $\lambda$ range over the finitely many eigenvalues of $A$, and similarly for $\mu, \nu$ we get all the simultaneous eigenspaces, and if you want an eigenbasis just choose any basis for each $V(\lambda, \mu, \nu)$ and take their union.