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$.
Here's a rough lower bound of the probability that $\Vert A\Vert<1$ given that $\Vert v\Vert=1$. Denote by $\sigma_k$ the area of the $k$-dimensional unit sphere.
Denote by $\lambda_\min$ respectively $\lambda_\max$ the smallest and respectively the largest of the eigenvalues of $A^2$. The nontrivial case is $\lambda_\min<1<\lambda_\max$, i.e., $\lambda_1=\lambda_\min$, $\lambda_n=\lambda_\max$.
Denote by $(v_1,\dotsc, v_n)$ the Euclidean coordinates determined by the basis $e_1,\dotsc, e_n$. Then $$ \Vert A v\Vert^2=\sum_{k=1}^n \lambda_kv_k^2. $$ if $\Vert v\Vert=1$, then $v_1^2=1-\sum_{k=2}^n v_k^2$ and $$ \Vert A v\Vert^2=\lambda_\min+\sum_{k=2}^n(\lambda_2-\lambda_\min)v_k^2 \leq\lambda_\min +(\lambda_\max-\lambda_\min)\sum_{k=2}v_k^2. $$ Thus, if $\Vert v\Vert =1$ and $$ \sum_{k=2}v_k^2<\frac{1-\lambda_\min}{\lambda_\max-\lambda_\min}=:\gamma(A), $$ then $\Vert Av\Vert<1$. In other words, if $\Vert v\Vert=1$ and $|v_1|>\sqrt{1-\gamma(A)}$, then $\Vert Av\Vert<1$. $\DeclareMathOperator{\Prob}{Prob}$. Hence $$ \Prob\big[ \Vert Av\Vert<1\big|\;\Vert v\Vert=1\big]\geq \frac{1}{\sigma_{n-1}}{\rm area}\;\big(\{v\in S^{n-1};\;v_1^2>1-\gamma(A)\}\big) $$ $$ =\frac{\sigma_{n-2}}{\sigma_{n-1}}\int_{1-\gamma(A)\leq v_1^2\leq 1}^1 (1-v_1^2)^{(n-3)/2} dv_1=\frac{\sigma_{n-2}}{\sigma_{n-1}}\int_{1-\gamma(A)}^1 (1-s)^{(n-3)/2}s^{-1/2} ds. $$ (The component $v_1$, as a random variable defined on the unit sphere has a Beta distribution) $$ \frac{\sigma_{n-2}}{\sigma_{n-1}}\int_{1-\gamma(A)}^1 (1-s)^{(n-3)/2}s^{-1/2} ds=\frac{1}{B((n-1)/2,1/2)}\int_0^{\gamma(A)} s^{(n-3)/2}(1-s)^{-1/2} ds. $$
We have $$ B((n-1)/2,1/2)\sim \Gamma(1/2)\sqrt{\frac{2}{n-1}}, $$ and $$ \int_0^{\gamma(A)} s^{(n-3)/2}(1-s)^{-1/2} ds \approx \frac{2}{n-1}\gamma(A)^{\frac{n-1}{2}}.$$
The expected time to detect $\Vert Av\Vert<1$ is about $\frac{1}{\sqrt{n}}\gamma(A)^{-n/2}$
Hence for $n<30$ the sample size ought to be in the millions. I cannot estimate how long it would take to run such samples. (I've run samples of this size on my laptop.) The computation of $\Vert Av\Vert$ is not very time consuming since $A$ is a $0/1$ matrix.
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.
Question. A natural question comes to mind. Suppose that $$ A=(A_{ij})_{1\leq i,j\leq n} $$ is a random symmetric $0/1$ matrix: the coefficients $A_{ij}$, $1\leq i\leq j\leq n$ are independent Bernoulli random variables with probability of success $1/2$. It would be interesting to know the probaility $p(n)$ that $A$ has an eigenvalue in $[-1,1]$
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.