Consider a determined linear system of differential order $k$ of the form $\sigma^{i_1\cdots i_k}_{ab}(x) \partial_{i_1} \cdots \partial_{i_k} u^b(x) + l.o.t = 0$. The coefficients $\sigma^{i_1\cdots i_k}_{ab}$, a square matrix in the $ab$ indices, constitute its _principal symbol_. By replacing $\partial_i$ by an indeterminate $p_i$ (like in the Fourier transform), the symbol defines a square matrix valued function on the cotangent space at $x$, $\sigma_{ab}(x,p) = \sigma^{i_1\cdots i_k}_{ab}(x) p_{i_1} \cdots p_{i_k}$, also referred to as the symbol. The locus $\mathcal{C}(x)$ of all $p$ for which the symbol fails to be an invertible matrix is the _characteristic variety_ of the PDE at $x$.
Further classification is based on the geometry of $\mathcal{C}(x)$, which is usually assumed to behave somewhat uniformly with respect to $x$. If $\mathcal{C}(x) = \{0\}$ then the equation is called _elliptic_. If $\mathcal{C}(x)$ looks like a cone which has compact intersection with some affine codimension-1 hyperplane (a "spacelike" hyperplane) then the equation is considered in a rather weak sense to be _hyperbolic_ (the geometric optics approximation could be used to construct waves traveling at finite speed). The simplest situation is when this intersection is smooth and is the boundary of some bounded convex set (e.g., a sphere). In general, the intersection can be multi-sheeted, non-convex, self-intersecting, etc. There is a gradation of notions of hyperbolicity corresponding to these variations.
In the usual case of a single scalar equation in two variables, the sign of the $A$-$B$-$C$-determinant can be used to determine the geometry of the characteristic variety. In principle, one could try to do something similar in higher dimensions, based on some more sophisticated invariants of the principal symbol coefficients $\sigma^{i_1\cdots i_k}_{ab}(x)$. But this will obviously become increasingly difficult in higher dimensions of dependent and independent variables.
Not much else can be said at this level of generality. But if you're interested in this level of generality, you'll find more information for instance in
> <cite authors="Seiler, Werner M.">_Seiler, Werner M._, [**Involution. The formal theory of differential equations and its applications in computer algebra**](http://dx.doi.org/10.1007/978-3-642-01287-7), Algorithms and Computation in Mathematics 24. Berlin: Springer (ISBN 978-3-642-01286-0/hbk; 978-3-642-01287-7/ebook). xxii, 650 p. (2010). [ZBL1205.35003](https://zbmath.org/?q=an:1205.35003).</cite>
----
**Update:** The following answers the updated question and the idea is to show how cumbersome it is to check the ellipticity condition for generic coefficients.
To check the ellipticity condition, you basically want to know whether the determinant $\tau(\xi,\eta)$ has the same sign everywhere outside $(\xi,\eta)\ne(0,0)$, lets say positive. You are hugely in luck! Because, according to a theorem of Hilbert (cf [Hilbert's 17th problem](https://en.wikipedia.org/wiki/Hilbert%27s_seventeenth_problem) and [positive semidefinite polynomials](https://en.wikipedia.org/wiki/Positive_polynomial)), the condition $\tau(\xi,\eta)\ge 0$ for a homogeneous polynomial in two variables is is one of the few cases when it is equivalent to the existence of a sum of squares representation $\tau(\xi,\eta) = \sum_i q_i(\xi,\eta)^2$. If a sum of squares representation of $\tau(\xi,\eta) = V\xi^4 + W\xi^3\eta + X\xi^2\eta^2 + Y\xi\eta^3 + Z\eta^4$ exists, then it must come from the diagonalization of a quadratic form.
Namely, there is a 1-parameter family of representations of $\tau(\xi,\eta)$ as a quadratic form
$$
\tau(\xi,\eta)
= \begin{pmatrix} \xi^2 & \xi\eta & \eta^2 \end{pmatrix}
\begin{pmatrix}
V & W & -\lambda \\
W & X+2\lambda & Y \\
-\lambda & Y & Z
\end{pmatrix}
\begin{pmatrix} \xi^2 \\ \xi\eta \\ \eta^2 \end{pmatrix} ,
$$
with $\lambda$ the only redundancy in this representation. So $\tau(\xi,\eta)$ is a sum of squares if and only if the above form is positive semidefinite for some value of $\lambda$. By [Sylvester's criterion](https://en.wikipedia.org/wiki/Sylvester%27s_criterion), the necessary and sufficient condition is the non-negativity of all of its principal minors:
$$\begin{aligned}
0 &\le V , \\
0 &\le X+2\lambda , \\
0 &\le Z , \\
0 &\le VX-W^2+2V\lambda , \\
0 &\le VZ - \lambda^2 , \\
0 &\le XZ-Y^2 +2Z\lambda , \\
0 &\le V X Z-W^2 Z - V Y^2 + 2 (V Z - 2 W Y) \lambda -X \lambda^2 -2 \lambda^3 .
\end{aligned}$$
If one wants, each of the above polynomial inequalities on $\lambda$ may be converted into an interval bound on $\lambda$, with the end points of the intervals explicitly expressed in terms of the coefficients of $\tau$, but the presence of quadratic and cubic polynomials will make these formulas rather cumbersome, involving several special cases depending on the discriminants of these polynomials.