I hope to use this answer to convince you that in general the method of characteristics cannot work for higher order PDEs in more than 2 variables. Nevertheless, there are some ideas in PDEs that are, in some sense, the next best thing, and maybe that will satisfy you.
First though: we will in this answer focus only on linear (systems) of differential equations, the philosophy being that if a tool works in general on nonlinear problems, then it should also work on linear problems.
Secondly, we can perform a reduction of order to focus only on first order systems of differential equations. A precise description of what you need to do is a bit annoying to write down, so I should just give an example instead: say you want to study the linear wave equation on $\mathbb{R} \times\mathbb{R}^2$, which we can write
$$ \partial^2_{tt} u - \partial^2_{xx} u - \partial^2_{yy} u = 0 $$
Then consider the vector-valued function
$$ U = (\partial_t u, \partial_x u, \partial_y u, u) $$
we see that it must solve a first order system
$$ \partial_t U = \begin{pmatrix} 0 & 1 & 0 & 0 \\
1 & 0 & 0 & 0 \\
0 & 0 & 0 & 0 \\
0 & 0 & 0 & 0\end{pmatrix} \partial_x U +
\begin{pmatrix} 0 & 0 & 1 & 0 \\
0 & 0 & 0 & 0 \\
1 & 0 & 0 & 0 \\
0 & 0 & 0 & 0\end{pmatrix} \partial_y U +
\begin{pmatrix} 0 & 0 & 0 & 0 \\
0 & 0 & 0 & 0 \\
0 & 0 & 0 & 0 \\
1 & 0 & 0 & 0\end{pmatrix} U $$
So, now, let us think in general about linear, constant coefficient partial differential systems of first order. Our domain is $\mathbb{R}^n$ (or a subset thereof). We shall assume $n \geq 2$ since otherwise we are already dealing with a system of ODEs. Our unknown is $U$ which takes values in $\mathbb{R}^k$. Since we have $k$ unknowns, we expect to have (for a well-determined problem) $k$ equations.
Our system can be written in the following form:
$$ 0 = \sum_{i = 1}^n \sum_{\beta = 1}^k A_{\alpha\beta}^i \partial_i U_\beta + \sum_{\beta = 1}^k B_{\alpha\beta} U_\beta $$
The method of characteristics, if it were to succeed, would reduce this system to a set of $k$ ordinary differential equations for components of $U$.
In other words, there would be a basis $v_1, \ldots, v_k$ of $\mathbb{R}^k$, such that, if we were to write $U = \sum u_\alpha v_\alpha$ in the basis, the coefficients would satisfy ODEs of the form
$$ \sum_{i = 1}^n \xi_i^{(\alpha)} \partial_i u_\alpha = \sum_{\beta = 1}^k \tilde{B}_{\alpha\beta} u_\beta $$
where for each $\alpha$, the object $\sum_{i} \xi_i^{(\alpha)} \partial_i$ is regarded as a vector field.
A little bit of linear algebra tells you then what this would require, is that the $n\times k \times k$ tensor $A^i_{\alpha\beta}$ can be simultaneously diagonalized. More precisely, we need that the $n$ matrices
$$ \mathbf{A}^{(i)}: (\mathbf{A}^{(i)})_{\alpha\beta} = A^i_{\alpha\beta} $$
can be simultaneously diagonalized, or that there exists a fixed set of left- and right-eigenvectors $r_1, \ldots, r_k$ and $\ell_1, \ldots, \ell_k$ such that each
$$ \mathbf{A}^{(i)} = \sum_{\beta = 1}^k \lambda_\beta^{(i)} \ell_\beta \otimes r_\beta. $$
Furthermore, for this to be representing the method of characteristics, the eigenvalues $\lambda_\beta^{(i)}$ must all be real, as they map to the coefficients $\xi_i^{(\alpha)}$ which needs to be interpreted as the characteristic vector fields on $\mathbb{R}^n$.
Here's the punch line For generic systems of differential equations, the matrices $\mathbf{A}^{(i)}$ are a set of $n$ generic matrices, and there is no hope that one can simultaneously diagonalize them.
However, in some special cases this diagonalization is possible.
- When the equation is scalar: in this case $k = 1$ and $1$-dimensional matrices are trivially diagonalized always. The very well developed method of characteristics for nonlinear first order scalar PDEs is what happens when you extend this as far as it takes you.
- When the the total number of dimensions is $n = 2$. In this case we just need to simultaneously diagonalize two matrices, which we get nice sufficient conditions from spectral theorems. For example, assuming that $\mathbf{A}^{(1)}$ is invertible, then a sufficient condition is that the matrix $(\mathbf{A}^{(1)})^{-1} \mathbf{A}^{(2)}$ is symmetric. This leads to the well-developed theory for "symmetric hyperbolic systems" in "1 space and 1 time dimension".
Note that in the second case, unlike the first, diagonalizability is not guaranteed. So for example, you cannot use the method of characteristics to study fundamentally elliptic problems. An example being the Cauchy-Riemann equations
$$ \partial_x U_1 = \partial_y U_2, \qquad \partial_x U_2 = - \partial_y U_1. $$
Trying to diagonalize this will lead to imaginary eigenvalues, which do not represent physical vector fields.
But all hope is not lost!
The fundamental idea of the method of characteristics is that it is easier to solve ODEs than PDEs. Given a PDE (especially a constant coefficient linear one), we can look for those solutions that are "essentially ODEs". In other words, we can pose the ansatz
$$ U(x) = \phi(\xi\cdot x) $$
where $\xi$ is some unit vector in $\mathbb{R}^n$ and $\phi:\mathbb{R}\to\mathbb{R}^k$ a function. These types of solutions will satisfy
$$ \sum_i \sum_{\beta} \xi_i A^i_{\alpha\beta} \phi'_\beta + \sum_\beta B_{\alpha\beta} \phi_\beta = 0 $$
which is a system of ODEs for $\phi$. Solutions of these types are called "plane waves".
A natural question then is: can we reconstitute a general solution $U$ as "linear combinations" of such plane wave solutions.
And this question is essentially what Fourier theory is all about.
From this point of view, the study of partial differential equations using harmonic analysis techniques is an extension of the method of characteristics. For linear but variable coefficients, this develops into the approach using "microlocal analysis" (and if you read the literature, a lot of the guiding principles of microlocal analysis is based on wave packets traveling along characteristic curves, and is very similar to what one would expect from the method of characteristics). If you are interested in these types of developments, the standard reference is Hormander's four volume "Analysis of Linear Partial Differential Operators".
