It's easier to understand if you change to new variables $\xi=x+ct$ and $\eta=x-ct$. Then the PDE becomes $\partial^2 u/\partial \xi \partial \eta = 0$. Integration with respect to $\eta$ gives $\partial u/\partial \xi = f(\xi)$ with $f$ an arbitrary function (the "constant" of integration, constant here meaning "independent of $\eta$"). Integrating again, now with respect to $\xi$, gives $u=F(\xi)+g(\eta)$ with $g$ an arbitrary function and $F$ the antiderivative of $f$.
$\newcommand{\+}{^{\dagger}}
\newcommand{\angles}[1]{\left\langle\, #1 \,\right\rangle}
\newcommand{\braces}[1]{\left\lbrace\, #1 \,\right\rbrace}
\newcommand{\bracks}[1]{\left\lbrack\, #1 \,\right\rbrack}
\newcommand{\ceil}[1]{\,\left\lceil\, #1 \,\right\rceil\,}
\newcommand{\dd}{{\rm d}}
\newcommand{\down}{\downarrow}
\newcommand{\ds}[1]{\displaystyle{#1}}
\newcommand{\expo}[1]{\,{\rm e}^{#1}\,}
\newcommand{\fermi}{\,{\rm f}}
\newcommand{\floor}[1]{\,\left\lfloor #1 \right\rfloor\,}
\newcommand{\half}{{1 \over 2}}
\newcommand{\ic}{{\rm i}}
\newcommand{\iff}{\Longleftrightarrow}
\newcommand{\imp}{\Longrightarrow}
\newcommand{\isdiv}{\,\left.\right\vert\,}
\newcommand{\ket}[1]{\left\vert #1\right\rangle}
\newcommand{\ol}[1]{\overline{#1}}
\newcommand{\pars}[1]{\left(\, #1 \,\right)}
\newcommand{\partiald}[3][]{\frac{\partial^{#1} #2}{\partial #3^{#1}}}
\newcommand{\pp}{{\cal P}}
\newcommand{\root}[2][]{\,\sqrt[#1]{\vphantom{\large A}\,#2\,}\,}
\newcommand{\sech}{\,{\rm sech}}
\newcommand{\sgn}{\,{\rm sgn}}
\newcommand{\totald}[3][]{\frac{{\rm d}^{#1} #2}{{\rm d} #3^{#1}}}
\newcommand{\ul}[1]{\underline{#1}}
\newcommand{\verts}[1]{\left\vert\, #1 \,\right\vert}
\newcommand{\wt}[1]{\widetilde{#1}}$
Lets $\ds{\xi \equiv x - ct\,,\quad \eta \equiv x + ct}$. Then,
\begin{align}
\partiald{}{x}&
=\partiald{\xi}{x}\,\partiald{}{\xi} + \partiald{\eta}{x}\,\partiald{}{\eta}
=\partiald{}{\xi} + \partiald{}{\eta}
\\[3mm]
\partiald{}{t}&
=\partiald{\xi}{t}\,\partiald{}{\xi} + \partiald{\eta}{t}\,\partiald{}{\eta}
=-c\,\partiald{}{\xi} + c\,\partiald{}{\eta}
\end{align}
\begin{align}
\partiald[2]{}{x}&=\partiald[2]{}{\xi} + 2\,{\partial^{2} \over \partial\xi\,\partial\eta} + \partiald[2]{}{\eta}
\\[3mm]
\partiald[2]{}{t}&=c^{2}\pars{%
\partiald[2]{}{\xi} - 2\,{\partial^{2} \over \partial\xi\,\partial\eta} + \partiald[2]{}{\eta}}
\end{align}
$$
\partiald[2]{}{x} - {1 \over c^{2}}\,\partiald[2]{}{t}
=4\,{\partial^{2} \over \partial\xi\,\partial\eta}
$$
$$
0=\partiald[2]{u}{x} - {1 \over c^{2}}\,\partiald[2]{u}{t}
=4\,{\partial^{2}u \over \partial\xi\,\partial\eta}
\quad\imp\quad\partiald{u}{\eta} = {\rm F}\pars{\eta}
\quad\imp\quad u = \overbrace{\int{\rm F}\pars{\eta}\,\dd\eta}^{\ds{\equiv\fermi\pars{\eta}}} + {\rm g}\pars{\xi}
$$
$$\color{#44f}{\large%
u = \fermi\pars{x + ct} + {\rm g}\pars{x - ct}}
$$
Best Answer
Maybe let's try first the 2D case: $u(x,y,t)$ and $u_{tt} = c^2\Delta (u_{xx} + u_{yy})$.
Then define $v^{(x)}_t = cu_x, v^{(y)}_t = cu_y$. For the spatial derivatives of $\boldsymbol v = \begin{pmatrix} v^{(x)} \\ v^{(y)} \end{pmatrix}$ require that $v^{(x)}_x + v^{(y)}_y = \frac 1c u_t$. Then for continuously differentiable $v^{(x)}, v^{(y)}$: \begin{align} \frac 1c u_{tt} &= \partial_t \Big(v^{(x)}_x + v^{(y)}_y \Big) \\ &= v^{(x)}_{xt} + v^{(y)}_{yt} \\ &= v^{(x)}_{tx} + v^{(y)}_{ty} \\ &= cu_{xx} + c u_{yy}\end{align} so that the original PDE (Wave equation) remains preserved. So the first order system reads $$ \partial_t\begin{pmatrix} v^{(x)} \\v^{(y)} \\ u \end{pmatrix} + \nabla \cdot \begin{pmatrix} -cu & 0 \\ 0 & -cu \\ -cv^{(x)} & -cv^{(y)} \end{pmatrix} = \begin{pmatrix} 0 \\ 0 \\ 0\end{pmatrix} $$ where the $\nabla \cdot$ acts row-wise.
This is the general form of a conservation law in multiple dimensions in divergence form: $$ \partial_t \boldsymbol u + \nabla \cdot \boldsymbol f(\boldsymbol u) = \boldsymbol{0}.$$
I am not sure how you would write a linear system in divergence form, since for a row-wise acting divergence $\nabla \cdot$ you need that $A \boldsymbol u \in \mathbb R^{m\times d}$, with $m$ the number of variables $| \boldsymbol u |$ and $d$ the spatial dimension. While you can ensure that $A$ has $m$ rows, there is no way for $A$ having $d$ columns if you multiply it with the column vector $\boldsymbol u \in \mathbb R^{m \times 1}$.
The extension to 3D is then straightforward, here you have $$ \partial_t\begin{pmatrix} v^{(x)} \\v^{(y)} \\ v^{(z)} \\ u \end{pmatrix} + \nabla \cdot \begin{pmatrix} -cu & 0 & 0\\ 0 & -cu &0 \\ 0 & 0 & -cu \\ -cv^{(x)} & -cv^{(y)} & -cv^{(z)} \end{pmatrix} = \begin{pmatrix} 0 \\ 0 \\ 0 \\ 0\end{pmatrix} $$