A symmetric matrix $A$ has an eigenvalue in $(-1,1)$ iff $A^2$ has an eigenvalue in $[0,1]$. Equivalently this means that $$ \min_{\Vert v\Vert=1} \frac{\Vert Av\Vert^2}{\Vert v\Vert^2} <1. $$ Hence you need to decide if the minimum of a nonnegative quadratic function on the unit sphere is $<1$.
I can think of two ways of doing this but I cannot comment on the computational complexity.
The first is a Monte Carlo type method. Sample $v$ uniformly on the sphere and compute $\Vert A v\Vert$. If the sample is large and there exists an eigenvalue of $A^2$ that is $<1$ you can pick detect this with high confidence. Technically, if the size of of the matrix is $n\times n$, then choose the coordinates $v_1,\dotsc, v_n$ of the random vector $v$ to be independent standard normal random variables and test if $$ \Vert A v\Vert^2 <\Vert v\Vert^2. $$
Here's a justifications. One can choose an orthonormal basis of your ambient space consisting of eigenvectors $e_1,\dotsc, e_n$ of $A^2$ corresponding to eigenvalues $\lambda_1\leq \cdots \leq \lambda_n$. If $\lambda_1<1$ then $\Vert Av\Vert^2<1$ in a neighborhood $U$ of $e_1$ on the unit sphere $\{\Vert v\Vert=1\}$. The probability that a random vector $v$ lands in $U$ is proportional to the size of $U$. One can provide lower bounds on the volume of $U$ (I won't do it here) that will give give you a lower bound on the probability of a random $v$ landing in $U$. That will give you and idea of the size of of the sample.
I cannot comment on the computational complexity of this approach.
If one your samples yield a vector such that $\Vert Av\Vert^2<1$ then game over.
If your samples yield only numbers $\Vert Av\Vert^2$ substantially bigger than $1$ then you can say with high confidence that $A$ has no small eigenvalues. If the numbers $\Vert A v\Vert$ are bigger than $1$ "many" are close to $1$ there is a bit of ambiguity.
The second approach is by gradient descent. Consider the function $$ f:\{\Vert v\Vert=1\}\to [0,\infty),\;\;f(v)=\Vert Av\Vert^2. $$ A flow line $v(t)$ of the negative gradient flow $$ \frac{dv}{dt}=-\nabla f(v) $$ will converge exponentially to an eigenvector of $A^2$. If the initial condition $v(0)$ is uniformly random on the unit sphere, then, with probability $1$, this flow line will converge exponentially to an eigenvector corresponding to a minimal eigenvalue.
Solve numerically this equation with random initial condition $v(0)$. I speculate that this discretization will lead you very fast to a decision concerning small e-values. The speed of convergence depends on how "packed" are the eigenvalues of $A^2$: if they are all packed in a small interval, the convergence will be slower.
Again, I cannot comment on the computational complexity of this approach.