Numerically solving a system of PDEs using change of variables. I am trying to numerically solve the system of PDE's written in matrix form as:
$$     \begin{bmatrix} (\lambda + 2\mu)\partial_x & \lambda\partial_y \\ \mu\partial_y & \mu\partial_x \end{bmatrix}\begin{pmatrix} u \\ v \end{pmatrix} = \begin{pmatrix} 0 \\ 0 \end{pmatrix}
 $$
using a finite difference discretization and multigrid Gauss-Seidel relaxation to solve the resulting linear system. In a course handout I read that this could be done efficiently by decoupling the equations using a change of variables using the transpose of the cofactor matrix of the system above:
$$ \begin{pmatrix} u \\ v \end{pmatrix} = \begin{bmatrix} \mu\partial_x & -\lambda\partial_y \\ -\mu\partial_y & (\lambda+2\mu)\partial_x \end{bmatrix} \begin{pmatrix} t_1 \\ t_2 \end{pmatrix}
 $$
so that the system decouples as:
$$ \begin{bmatrix} (\lambda\mu + 2\mu^2)\partial_x^2 - \lambda\mu\partial_y^2 & 0 \\ 0 & (\lambda\mu + 2\mu^2)\partial_x^2 - \lambda\mu\partial_y^2 \end{bmatrix} \begin{pmatrix} t_1 \\ t_2 \end{pmatrix} = \begin{pmatrix} 0 \\ 0 \end{pmatrix} $$
Now I get that I could discretize this, solve for $t_1$ and $t_2$ and then compute $u$ and $v$. However in the handout it is stated that a full change of variables is not actually necessary, and that "the solution process can be arranged such" that I am effectively solving the decoupled system, but in terms of the original variables $u$ and $v$.
Unfortunately the explanation ends there. How would this work in practice? 
 A: Rewrite the system as $AU_x + BU_y =0$,
where $U =(u,v)^T$ and
$$
A=\begin{bmatrix}
\lambda +2\mu &0\\
0& \mu
\end{bmatrix} ,\qquad
B=\begin{bmatrix}
0 &\lambda\\
\mu& 0
\end{bmatrix} .
$$
For non-zero parameters, $A$ is invertible and the system reads
$U_x + M U_y = 0$
with $M=A^{-1}B$. This system can be solved numerically and analytically by using diagonalization, and then transforming back to $U$. Setting $V=P^{-1}U$ and $M = P\Lambda P^{-1}$ with
$$
P=\begin{bmatrix}
-c & c\\
1 & 1
\end{bmatrix} ,\qquad
P^{-1}=\tfrac12\begin{bmatrix}
-1/c & 1\\
1/c & 1
\end{bmatrix}, \qquad
c = \sqrt{\tfrac{\lambda}{\lambda+2\mu}}\; ,
$$
the system is rewritten as $V_x + \Lambda V_y = 0$ where $\Lambda = \text{diag}\lbrace{-c}, c\rbrace$ is diagonal.
A standard numerical procedure consists in using finite-volume schemes, such as the Upwind scheme, the Lax-Friedrichs method or the Lax-Wendroff method.
The analytical procedure is based on the method of characteristics (several related posts on this site), which provides the coordinates $(F(y+cx), G(y-cx))^T$ of $V$ where $F$, $G$ are arbitrary functions deduced from the boundary conditions. The coordinates of $U = PV$ are therefore
\begin{aligned}
u(x,y) &= c \big( G(y-cx) - F(y+cx)\big) \\
v(x,y) &= G(y-cx) + F(y+cx)
\end{aligned} 
