What follows below is a brief description of the procedure Dirac used to obtain the Dirac equation. For my question, which is about this procedure, see the last paragraph.
Historically, Dirac was trying to find a relativistic version of the Schrödinger equation, just like the Klein–Gordon equation,
$$ \tag{*}\label{*} \left(-\frac{1}{c^2}\frac{\partial^2}{\partial t^2} + \nabla^2\right)\phi = \frac{m^2c^2}{\hbar^2}\phi, $$
but which was of first order in space and time instead of second order. (This is because while the Klein–Gordon equation is Lorentz invariant and describes a particle that behaves like a particle described by the Schrödinger equation in low velocities, it is too weak in the sense that it also allows solutions with non-conserved probability currents if the normal probability density $\rho = \phi^*\phi$ is used, or alternatively, if one changes to a probability density current that is conserved, solutions with negative probability densities.)
He proceeded to take the square root of the wave operator (the parenthesis in the LHS of $\eqref{*}$), and assert that this was of firs order, thus:
$$ \tag{**}\label{**} \nabla^2 - \frac{1}{c^2}\frac{\partial^2}{\partial t^2} = \left(A \partial_x + B \partial_y + C \partial_z + \frac{i}{c}D \partial_t\right)\left(A \partial_x + B \partial_y + C \partial_z + \frac{i}{c}D \partial_t\right). $$
The right hand side of this equation expands to
$$ \begin{array}{c} \displaystyle A^2 \partial_x^2 + B^2 \partial_y^2 + C^2 \partial_z^2 + \frac{-1}{c^2}D^2 \partial_t^2 \\ \displaystyle +\, (AB+BA)\partial_x\partial_y + (AC+CA)\partial_x\partial_z + \left(A\frac{i}{c}D+\frac{i}{c}DA\right)\partial_x\partial_t + (BC+BA)\partial_y\partial_z + \left(B\frac{i}{c}D+\frac{i}{c}DB\right)\partial_y\partial_t + \left(C\frac{i}{c}D+\frac{i}{c}DC\right)\partial_z\partial_t. \end{array} $$
For this to match the LHS of $\eqref{**}$, we require that $A^2=B^2=C^2=D^2=1$. However, at the same time we require that $AB+BA=AC+CA=AD+DA=BC+CB=BD+DB=CD+DC=0$, which if $A,B,C,D$ are real or complex numbers is impossible, since the requirements imply that $A$, $B$, $C$ and $D$ must not commute. However, this is possible if $A$, $B$, $C$ and $D$ are matrices. We already know that satisfying these requirements is possible for three matrices $A$, $B$ and $C$ (these are the Pauli matrices, which are three $2\times 2$ matrices). For four matrices, it turns out that these need to be of at least size $4\times 4$, for example the gamma (or Dirac) matrices $\gamma^1$, $\gamma^2$, $\gamma^3$ and $\gamma^0$ multiplied by $i$. Hence, by taking the square root of the factor before $\phi$ on both sides of $\eqref{*}$, we get
$$ i\left(\gamma^1 \partial_x + \gamma^2 \partial_y + \gamma^3 \partial_z - \frac{1}{c}\gamma^0 \partial_t\right)\phi = \frac{mc}{\hbar}\phi, $$
which is known as the Dirac equation. (By using the convention $c=\hbar=1$, as well as the Feynman slash notation, this can also be written as $\left(i\partial\!\!\!/ \ \,- m\right)\phi = 0$.) However, as we see, $\phi$ is now required to bee a vector-valued function (more specifically, $\phi: \mathbb{R}^d\to \mathbb{C}^4$, since the gamma matrices are of size $4\times 4$, where $d=4$ is the dimensionality of the spacetime) instead of a complex-valued function, and even if $\phi$ would be vector-valued in $\eqref{*}$ as well, that equation would be weaker than the Dirac equation since it would allow solutions which the Dirac equation doesn't allow.
My question is, is this procedure to reduce the order of a partial differential equation known by any specific name? And has it been used in other contexts as well, except for turning the Klein–Gordon equation into the Dirac equation?