The numerical range of a real matrix is an interval

Step 1: only the symmetric part matters. $x^TAx$ is a real number, so it equals its transpose $x^TA^Tx$. Hence $x^TAx=\frac12(x^TAx+x^TA^Tx)=x^TSx$.

Step 2: the endpoints. $S$ is symmetric; let $\lambda_1\le\dots\le\lambda_n$ be its eigenvalues with an orthonormal eigenbasis $v_1,\dots,v_n$. For a unit vector $x=\sum c_iv_i$ we have $\sum c_i^2=1$ and
$$x^TSx=\sum_i\lambda_ic_i^2\in[\lambda_1,\lambda_n],$$
with equality at $x=v_1$ and $x=v_n$. So $\min$ and $\max$ of $x^TSx$ on the sphere are $\lambda_{\min}$ and $\lambda_{\max}$.

Step 3: no gaps. The function $x\mapsto x^TSx$ is continuous. For $n\ge2$ the unit sphere $S^{n-1}$ is connected, so its image is a connected subset of $\mathbb{R}$, i.e. an interval, and it is compact. An interval containing $\lambda_{\min}$ and $\lambda_{\max}$ and contained in $[\lambda_{\min},\lambda_{\max}]$ is the whole of $[\lambda_{\min},\lambda_{\max}]$. (For $n=1$ the sphere is $\{\pm1\}$ and both sides equal $\{a\}$.)

Explicit path. For a direct proof, $x(t)=\cos t\,v_1+\sin t\,v_n$ gives $x(t)^TSx(t)=\lambda_1\cos^2t+\lambda_n\sin^2t$, which sweeps $[\lambda_1,\lambda_n]$ as $t\in[0,\pi/2]$.

Remark. For complex matrices the numerical range is a convex set in $\mathbb{C}$ (Toeplitz-Hausdorff theorem), not an interval. The real quadratic form only sees $S$, which is why a skew-symmetric $A$ has numerical range $\{0\}$.