Numerical Solution of Systems of PDE Could someone give me some reference  to the Numerical solution of a System of PDEs  of following type.. (which also encompasses strongly elliptic system of PDEs) in 2D or 3D.
$$\left\{
\begin{aligned}
& \mathrm{div} (\mathbf{U}_k)=f_k ,  \text{ for } k=1,\ldots,M 
\\ 
&\mathbf{U}_k = \sum^M_{l=1}B_{kl} \cdot \mathrm{grad\,}  P_l ,\text{ for } k=1,\ldots,M
\end{aligned}\right.
$$
where $B_{kl}$ is a matrix of appropriate dimension.
 A: $\newcommand{\dualp}[2]{{\langle {#1},{#2} \rangle}}$
$\newcommand{\v}[1]{\boldsymbol{#1}}$
As far as I know, if you are using finite element methods on an unstructured mesh, there are two types of approach available, suppose we have a polygonal/polyhedral domain $\Omega$ triangulated by a set of triangles/tetrahedra in 2D/3D respectively, and proper boundary conditions are assumed for the PDEs.


*

*First is mixed finite element method(MFEM), the numerical method for strongly elliptic equation system shares the same core idea with the MFEM for diffusion equation $-\mathrm{div}(D\, \mathbf{grad} \,u) = f$, rewrite this second order PDE as a first order system:
$$\left\{
\begin{array}{l}
&D^{-1}\,\boldsymbol{\sigma} + \mathbf{grad} \,u &= 0 
\\
& \mathrm{div} \boldsymbol{\sigma} &=f
\end{array}\right.
$$
which falls into the framework of the following saddle point problem:
$$\left\{
\begin{array}{l}
&A\boldsymbol{\sigma} + B^T \,u &= \v{g} &\text{ in } V'
\\
& B\boldsymbol{\sigma} &=f &\text{ in } Q'
\end{array}\right.
$$
where $V$ and $Q$ are two Hilbert spaces, in this case $V = H(\mathrm{div},\Omega)$ and $Q = L^2(\Omega)$. The variational formulation is 


$$\left\{
\begin{array}{l}
&\dualp{A\v{\sigma}}{\v{\tau}} + \dualp{\v{\tau}}{B^T \,u} &= \dualp{\v{g}}{\v{\tau}} &\forall \v{\tau}\in V
\\
& \dualp{B\boldsymbol{\sigma}}{q} &=\dualp{f}{q} &\forall q\in Q
\end{array}\right.
$$
Take the test function space be a finite dimensional space $V_h$, either Raviart-Thomas or Brezzi-Douglas-Marini element(conforming) or other nonconforming spaces, and you have the MFEM formulation, for detail please refer to Mixed and Hybrid Finite Element Methods by Brezzi and Fortin, in your case it is:
$$\left\{
\begin{array}{l}
&\dualp{B^{-1}\mathbf{U}_k}{\mathbf{V}_k} - \dualp{\mathbf{V}_k}{\mathbf{grad} P} &= 0 &\forall \,\mathbf{V}_k
\\
& \dualp{\mathrm{div}{\mathbf{U}_k}}{q_k} &=\dualp{f_k}{q_k} &\forall q_k
\end{array}\right.
$$


*

*Second is least-square finite element method(LSFEM), Google scholar gives me this as the most cited paper: First-Order System Least Squares for Second-Order Partial Differential Equations, and a recent review is this. The idea is, if we still use the diffusion equation as example, we would like to minimize the following least-square functional over appropriate finite element spaces:
$$
\mathcal{J}(v,\v{\tau}) = \|\v{\tau}+D\mathbf{grad}v\|^2_{\boldsymbol{L}^2(\Omega)} + \|\mathrm{div}\v{\tau} - f\|^2_{L^2(\Omega)}
$$
In your case, the LS functional you would like to minimize is 
$$
\mathcal{J}(\mathbf{V}_1,\ldots,\mathbf{V}_M,Q_1,\ldots,Q_M) = \sum_k \|\mathbf{V}_k - D\mathbf{grad}Q_k\|^2_{\boldsymbol{L}^2(\Omega)} + \sum_k \|\mathrm{div}\mathbf{V}_k - f_k\|^2_{L^2(\Omega)}
$$

