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.