You can apply a common approach known as the spectrum slicing. Assume that $M$ is a Hermitian matrix factored as
$$\tag{1}M=LDL^*$$ (a complex variant of the $LDL^T$ factorisation), where $L$ is unit lower triangular and $D$ diagonal. The inertia (the number of negative, zero, and positive eigenvalues) of $M$ and $D$ is the same and since $D$ is diagonal, you can get the inertia of $M$ pretty easily.
Now fix a $\sigma\in\mathbb{R}$. The eigenvalues of $M-\sigma I$ (where $I$ is the identity matrix) are the eigenvalues of $M$ minus $\sigma$. Consider the factorisation
$$
M-\sigma I=L_{\sigma}D_{\sigma}L_{\sigma}^*.
$$
Both $M-\sigma I$ and $D_{\sigma}$ have again the same inertia so, e.g., the number of negative diagonal entries of $D_{\sigma}$ determines the number of eigenvalues of $M$ smaller than $\sigma$.
In your case, since you know that the spectrum of $M$ is contained in the interval $[-1,1]$, you need to perform only two factorisations to determine how many eigenvalues are contained in the intervals $[-1,-\alpha)$ and $(\alpha,1]$ for some $\alpha\in(0,1)$. In particular,
$$
M+\alpha I=L_+D_+L_+^*, \quad\text{and}\quad M-\alpha I=L_-D_-L_-^*.
$$
Then, the number of negative diagonal entries in $D_+$ gives you the number of eigenvalues of $M$ smaller than $-\alpha$ {that is, in the interval $[-1,-\alpha)$} and the number of positive diagonal entries of $D_-$ gives the number of eigenvalues of $M$ greater than $\alpha$ {that is, in the interval $(\alpha,1]$}.
In MATLAB, the factorisation you look for is implemented in the function ldl.
Actually, you could also use eigs to compute, say, $k$ eigenvalues largest in magnitude (that is, close to the boundaries of the interval $[-1,1]$) and check if the result gave you an eigenvalue with a magnitude smaller than $\alpha$ with both signs). If it is so, then the result of eigs contains all eigenvalues in $[-1,-\alpha)$ and $(\alpha,1]$ plus some eigenvalues outside of these intervals. If this criterion is not satisfies, you need to increase $k$. Of course, this could be rather costly especially if the eigenvalues would be clustered around $-1$ and $1$ which would require to use relatively high $k$.
P.S.: The matrix $D$ is not generally diagonal, but block diagonal with $1\times 1$ and $2\times 2$ blocks. Still, the eigenvalues of $D$ (and hence their signs which matter here) can be easily obtained.
Best Answer
The difference between a full and a sparse matrix in MATLAB is not really a mathematical one (because all calculations should work correctly either way), but more a matter of how the entries of the matrix are stored.
The MATLAB documentation links to this paper for more information and the following example is given there:
The command
generates a full matrix, which can be displayed as
This means that all the zeros are actually stored explicitly in memory, because also the number $0$ needs the same amount of memory as every other number does.
If the full matrix is transformed to a sparse matrix with the command
it is displayed (and stored in memory) as
You can see here that only the nonzero entries of the matrix and in addition their positions inside the matrix are stored, whereas the zero elements are ignored and all entries of the matrix which are not present in the $(i,j)$-coordinates are implicitly assumed to be zero.
If a matrix has a large number of zeros, this way of storing it saves a lot of memory: If you have a $n\times n$ -matrix with $m$ nonzero elements, it needs $n^2$ units of memory to be stored as a full matrix, but $3m$ units of memory to be stored as a sparse matrix (assuming that storing the positions $(i,j)$ needs 2 units of memory). If $m$ is much smaller than $n^2$, i.e. the matrix is sparse, then you will see that $3m$ is of course also much smaller than $n^2$.