Solution (source code)

= Solution

Each coordinate of the transformed vectors evolves independently because $D$ and $E$ are <diagonal matrices>. Reorder the unknowns by collecting the $k$th components across all grid columns:
$$
\widehat v_k=((v_1)_k,\ldots,(v_m)_k)^T,
\qquad
\widehat c_k=((c_1)_k,\ldots,(c_m)_k)^T.
$$
Then system (3) splits into the $m$ uncoupled systems
$$
\boxed{\Lambda_k\widehat v_k=\widehat c_k,
\qquad k=1,\ldots,m,}
$$
where
$$
\boxed{\Lambda_k=[e_k,d_k,e_k]
=\begin{pmatrix}
d_k&e_k&&\\
e_k&d_k&\ddots&\\
&\ddots&\ddots&e_k\\
&&e_k&d_k
\end{pmatrix}.}
$$
\b[Thus the original $m^2\times m^2$ problem reduces to one <discrete sine transform> followed by $m$ independent tridiagonal solves and an inverse sine transform.]