The first thing to note is that the correspondence between finding
roots of a polynomial (any polynomial) and finding the eigenvalues
of an arbitrary matrix is really direct, and it's a rich subject, see Pseudozeros of
polynomials and pseudospectra of companion matrices by Toh and
Trefethen and the references
there.
Basically, the 2×2 case is trivial and the standard formula,
$$ x_1 = \frac{-b-\mathrm{sign}(b)\sqrt{Δ}}{2a}, \qquad x_2 = c/(ax_1), \qquad Δ = \det \begin{pmatrix} b & 2a\\ 2c & b \end{pmatrix} $$
is numerically stable and accurate, so long as the determinant
$Δ$ is evaluated accurately—the direct formula will be inaccurate when $b^2≈4ac$, but there is an accurate formula due to Kahan that uses FMAs (https://hal.inria.fr/ensl-00649347v1/document).
There is no such direct equivalent even for cubic polynomials. (Edit: see the link to Kahan's method below in CADJunkie's comment. This might well be wrong.) The
direct formula is not always numerically stable (contrary to your assumption, I think), and it can't be made
numerically stable the same way as the quadratic formula by inserting
the right signs somewhere. You could try evaluating it with extra precision, e.g. with double-native arithmetic. But the approaches that do work directly on
the polynomial are pretty complicated. For example
(https://doi.org/10.1145/2699468, which also works on quartic polynomials) you could use Newton's method with a
good precomputed first guess, but it gets really rather complicated, and the
speedup is not even all that large.
The explicit formulas for degree-4 polynomials are likewise not always
numerically stable. The polynomials that are the hardest tend to have
uncommon roots (small, or close to each other, differing in magnitude
by a lot), but even testing your code on a few billion purely random polynomials can usually reveal numerical errors.
One curious thing related to this is that Jenkins-Traub, which is a good common way of finding polynomial roots, is actually an eigenvalue algorithm (inverse iteration) in disguise.
I would say the explicitness of the formulas is, in a way, misleading:
it tricks you into thinking that because the formula has a closed form
that means it's somehow cheaper/easier. I would really recommend that you actually test/benchmark this on some test data.
It doesn't have to be true:
determining the roots of degree-$≥3$ polynomials is within a small
integer factor of difficulty of the full problem of determining the
eigenvalues of a small matrix, and the standard library routines for
eigenvalues are much more robust and well-tested. So by reducing the small
eigenvalue problem to a low-degree polynomial roots problem, you
aren't necessarily simplifying it.
What are the pros and cons to using the characteristic polynomial for getting eigenvalues this specifically for this case?
I think the main con is that this assumption that you make:
On the other hand, the equation is only quartic at best and we have analytical formulas for the polynomial roots so we shouldn't get too far off.
that because a formula is in closed form, and analytical, that means it's easy/cheap/accurate, is not necessarily true. It could be true on specific data you might have, but as far as I know it's not true in general.
P.S. The whole distinction between in-closed-form and not-in-closed-form gets really finicky with computer arithmetic: you might think that $\cos (\cdots)$ in a cubic formula is a closed form, but as far as computer arithmetic is concerned, that's just another approximate rational function, it is maybe faster, but not fundamentally different, from the approximate rational function that defines the result of an eigenvalue algorithm.