Let $W = M + D$, where $M$ is the original $n \times n$ matrix and $D$ is the added diagonal matrix that we want to determine.
$W$ is symmetric, thus diagonalizable by an adjoint action of the orthogonal group.
Larger multiplicities in the eigenvalues of $W$ imply smaller dimensions of the adjoint orbits.
For example, if we have an eigenvalue of multiplicity $n-1$, then the adjoint orbit will be
$O(n)/(O(n-1) \times O(1)) = S^{n-1}/Z_2 = RP^{n-1}$ of dimension $n-1$;
while, if all the eigenvalues are distinct, the adjoint orbit is the real flag manifold $Fl_{\mathbb{R}}^n = O(n)/(Z_2)^n$ of dimension $ \frac{n(n-1)}{2}$. (Of course, in the case of scalar multiple of the unit matrix, the adjoint orbit is just a single point).
Thus in order to achieve large multiplicities, we need to minimize the dimension of the adjoint orbit.
This problem can be reduced to a problem of matrix rank minimization as follows:
Let $\{l_i \}_{i=1,...,n(n-1)/2}$, be a set of generators of the Lie algebra of $O(n)$ normalized according to:
$\textrm{tr}(l_i l_j) = \delta_{ij}$. The dimension of the adjoint orbit equals the rank of the Gram matrix $C$ whose elements are given by, $C_{ij} = \textrm{tr}([l_i, W][l_j, W])$. The problem is thus reduced to the minimization of the rank of the Gram matrix whose elements depend linearly on the added diagonal matrix elements.
One of the possible methods to solve this problem is through a convex programming
heuristic approach for the solution of matrix rank minimization,
based on replacing the rank by the nuclear norm (the sum of the singular values), as explained in the following lecture notes by: P.A. Parillo. The nuclear norm is a convex envelope of the rank which may explain why this method works well in practice in general.