This is quite an interesting question, perhaps a research problem.
I think an elementary answer should be a high school algebra answer in the
sense of Abhyankar and it would have to be in the spirit of what follows.
But first a little story.

I was teaching linear algebra and had just covered eigenvalues and characteristic polynomials but was not yet at the chapter on the spectral theorem for real symmetric matrices. I was looking for problems to assign for my students as homework in the textbook we were using.
One of the exercises was to show that a real matrix
$$
A=\left[
\begin{array}{cc}
\alpha & \beta \\\
\beta & \gamma
\end{array}
\right]
$$
only had real eigenvalues.
Not too hard. Write
the characteristic polynomial
$$
\chi(\lambda)=det(\lambda I-A)=\lambda^2-(\alpha+\gamma)\lambda+\alpha\gamma-\beta^2
$$
then its discriminant is
$$
\Delta=(\alpha+\gamma)^2-4(\alpha\gamma-\beta^2)=(\alpha+\gamma)^2+4\beta^2\ge 0\ .
$$
Hence two real roots.

The next problem in the book was to do the same for
$$
A=\left[
\begin{array}{ccc}
\alpha & \beta & \gamma\\\
\beta & \delta & \varepsilon \\\
\gamma & \varepsilon & \zeta
\end{array}
\right]
$$
and (silly me) I also assigned it...

Here is the solution in the 3X3 case. All roots are real if the discriminant (for a binary cubic) is nonnegative. The discriminant of the characteristic polynomial is
$$
\Delta  = (\delta \varepsilon ^{2} + \delta \zeta ^{2} - 
\zeta \delta ^{2} - \zeta \varepsilon ^{2} + \zeta \alpha 
^{2} + \zeta \gamma ^{2} - \alpha \gamma ^{2} - \alpha 
\zeta ^{2} + \alpha \beta ^{2} + \alpha \delta ^{2} - \delta 
\alpha ^{2} - \delta \beta ^{2})^{2} \\\
\mbox{} + 14(\delta \gamma \varepsilon  - \beta 
\varepsilon ^{2} + \beta \gamma ^{2} - \alpha \gamma 
\varepsilon )^{2} \\\
\mbox{} + 2(\delta \alpha \gamma  + \delta \beta 
\varepsilon  + \delta \gamma \zeta  - \gamma \delta ^{2} - 
\gamma \varepsilon ^{2} + \gamma ^{3} - \alpha \beta 
\varepsilon  - \alpha \gamma \zeta )^{2} \\\
\mbox{} + 2(\delta \beta \gamma  + \delta \varepsilon 
\zeta  - \varepsilon ^{3} + \varepsilon \alpha ^{2} + 
\varepsilon \gamma ^{2} - \alpha \beta \gamma  - \alpha 
\delta \varepsilon  - \alpha \varepsilon \zeta )^{2} \\\
\mbox{} + 2(\zeta \alpha \beta  + \zeta \beta \delta 
 + \zeta \gamma \varepsilon  - \beta \varepsilon ^{2} - 
\beta \zeta ^{2} + \beta ^{3} - \delta \alpha \beta  - 
\alpha \gamma \varepsilon )^{2} \\\
\mbox{} + 14(\zeta \beta \varepsilon  - \gamma 
\varepsilon ^{2} + \gamma \beta ^{2} - \alpha \beta 
\varepsilon )^{2} \\\
\mbox{} + 2(\zeta \beta \gamma  + \delta \varepsilon 
\zeta  - \varepsilon ^{3} + \varepsilon \alpha ^{2} + 
\varepsilon \beta ^{2} - \alpha \beta \gamma  - \alpha 
\delta \varepsilon  - \alpha \varepsilon \zeta )^{2} \\\
\mbox{} + 14(\varepsilon \beta ^{2} + \zeta \beta \gamma 
 - \delta \beta \gamma  - \varepsilon \gamma ^{2})^{2} \\\
\mbox{} + 2(\zeta \alpha \beta  + \zeta \beta \delta 
 + \zeta \gamma \varepsilon  - \beta \gamma ^{2} - \beta 
\zeta ^{2} + \beta ^{3} - \delta \alpha \beta  - \delta 
\gamma \varepsilon )^{2} \\\
\mbox{} + 2(\alpha \gamma \zeta  + \zeta \beta 
\varepsilon  - \gamma ^{3} + \gamma \beta ^{2} + \gamma 
\delta ^{2} - \delta \alpha \gamma  - \delta \beta 
\varepsilon  - \delta \gamma \zeta )^{2}\ .
$$

This formula comes from a paper by Ilyushechkin in Mat. Zametki, **51**, 16-23, 1992.


I suspect the elementary answer should be as follows.
First find a list of invariants or covariants of binary forms $C_1,C_2,\ldots$
such that a form with real coefficients has only real roots iff these covariants are nonnegative. Apply this to the characteristic polynomial of a general real syymetric matrix and show you get sums of squares. I suppose these covariants, via Sturm's sequence type arguments, should correspond to subresultants or rather subdiscriminants.
This seems also related to Part 2) of Godsil's answer.