Constructing a Distributional Solution to the Inhomogeneous C.R. Equations The question is to find a fundamental solution to the system of equations in $\mathbb{R}^{2}$
\begin{array}{l}
u_{x}-v_{y}=f\\
u_{y}+v_{x}=g\end{array}
and to express the answer as a $2\times2$ matrix of tempered distributions.  I really have had no experience with solving systems of PDE, so while the Fourier transform will undoubtedly be needed here, I'm not sure how to apply it.  Also, I don't really understand the directions.  How can we have a $2\times2$ matrix of solutions? Are we not solving for $u$ and $v$?  Furthermore, what characterizes a fundamental solution for a system.  It almost certainly should involve the $\delta$ distribution, no?  So if this is the case, why the $f$ and $g$?

Here's my attempted solution...
I will use the following convention for the Fourier transform: $$\mathcal{F}(u(x,y))=\frac{1}{2\pi}\int\int_{\mathbb{R}^{2}}u(x,y)e^{-i(x\xi+y\zeta)}\;dxdy.$$  Then applying $\mathcal{F}$ to both equations we get
\begin{array}{l}
i\xi\hat{u}-i\zeta\hat{v}=\hat{f}\\
i\zeta\hat{u}+i\xi\hat{v}=\hat{g}.
\end{array}
If we solve for $\hat{u}$ and $\hat{v}$ we then get
\begin{array}{l}
\hat{u}=-\frac{i\xi\hat{f}+i\zeta\hat{g}}{\xi^{2}+\zeta^{2}}=-\frac{\hat{f_{x}}+\hat{g_{y}}}{\xi^{2}+\zeta^{2}}\\
\hat{v}=-\frac{i\zeta\hat{f}-i\xi\hat{g}}{\xi^{2}+\zeta^{2}}=-\frac{\hat{f_{y}}-\hat{g_{x}}}{\xi^{2}+\zeta^{2}}\end{array}
so that by the convolution theorem
\begin{array}{l}
u=-\frac{1}{2\pi}(f_{x}+g_{y})*\mathcal{F}^{-1}\left(\frac{1}{\xi^{2}+\zeta^{2}}\right)\\
u=-\frac{1}{2\pi}(f_{y}-g_{x})*\mathcal{F}^{-1}\left(\frac{1}{\xi^{2}+\zeta^{2}}\right)
\end{array}
So then, formally speaking, if
$$h=\frac{-1}{2\pi}\mathcal{F}^{-1}\left(\frac{1}{\xi^{2}+\zeta^{2}}\right),$$
then
$$u=f_{x}*h+g_{y}*h$$
and
$$v=f_{y}*h-g_{x}*h.$$

The fundamental solution is convolution with $f$ and $g$, however, not the partial derivatives of $f$ and $g$.  Also, the above solution can't be expressed with respect to convolution of a $2\times2$ matrix.  So, let us use the multipliers for the differentiation operators and immediately apply the convolution theorem.  To that end we have
$$\begin{array}{l}
\hat{u}=\frac{-i}{\xi^{2}+\zeta^{2}}\left(\xi\hat{f}+\zeta\hat{g}\right).\\
\hat{v}=\frac{-i}{\xi^{2}+\zeta^{2}}\left(\zeta\hat{f}-\xi\hat{g}\right).
\end{array}$$
Then if $$A=\frac{1}{2\pi i}\mathcal{F}^{-1}\left(\frac{\xi}{\xi^{2}+\zeta^{2}}\right),$$ $$B=\frac{1}{2\pi i}\mathcal{F}^{-1}\left(\frac{\zeta}{\xi^{2}+\zeta^{2}}\right),$$
and $$M=\begin{pmatrix}A&B\\B&-A\end{pmatrix}$$
then
$$\begin{pmatrix}u\\v\end{pmatrix}=M*\begin{pmatrix}f\\g\end{pmatrix}.$$
Would you think this is an acceptable answer?  Is it possible to explicitly invert $A$ and $B$?  The instructions ask for a $2\times2$ matrix of tempered distributions representing the fundamental solution. Clearly $M$ is the fundamental solution.  Presumably $f$ and $g$ are tempered so the operations used to solve for $u$ and $v$ are also tempered.  I still to this day just have difficulty accepting results obtained by formal operations on distributions as if they were functions as being true (or as being obtained rigorously I should say) since behind the scenes they are not functions (of course, I could probably just write everything as a linear functional to verify what I got already, but this seems dumb and unnecessary).
 A: In complex notation, you are asked to find a function $\varphi = u+iv$ such that $\frac{\partial \varphi}{\partial z} = \frac12(f+ig)$. The relevant formula from complex analysis (Cauchy-Pompeiu) is usually stated in terms of $\frac{\partial \varphi}{\partial z}$ derivative. Rather than tweak the formula, you can look for $\psi = u-iv$, which satisfies  $\frac{\partial \psi}{\partial \bar z} = \frac12(f+ig)$. A solution (non-unique) is given by 
$$
(u-iv)(z) = -\frac{1}{2}\iint \frac{f(\zeta)+ig(\zeta)}{z-\zeta} \,d\lambda
\tag{1}$$
where $\lambda$ is the 2-dimensional Lebesgue measure. To turn (1) into a real-variable formula, you'll have to separate the real and imaginary part. It's easier to work with convolution notation:
$$
(f+ig)*\frac{1}{z} = f*\frac{x}{x^2+y^2} + g* \frac{y}{x^2+y^2} +i\left( -f*\frac{y}{x^2+y^2} +g*\frac{x}{x^2+y^2} \right)$$
This is how you turn (1) into 
$$\begin{pmatrix} u \\ v\end{pmatrix} = A*\begin{pmatrix} f \\ g\end{pmatrix}$$ with some $2\times 2$ matrix $A$.
A: (In this post $\mathbb{K}$ stands for either the real or the complex field). 
Regarding the "fundamental solution" concept and the $\delta$ distribution. The linear operator that we are dealing with here is the following: 
$$
L \left(\begin{bmatrix} u \\v  \end{bmatrix}\right) = \begin{bmatrix} u_x-v_y \\u_x+v_y \end{bmatrix},
$$
and is defined on $\left[\mathscr{D}'(\mathbb{R}^2)\right]^2$. We can introduce a convolution product on this space via the following formula: 
$$\tag{1}
\begin{bmatrix} 
u_1 \\ v_1\end{bmatrix} \star 
\begin{bmatrix} 
u_2 \\ v_2 \end{bmatrix} = \begin{bmatrix}
u_1\ast u_2 + v_1\ast v_2 \\
u_1 \ast v_2 + v_1\ast u_2
\end{bmatrix}^{[1]},
$$
where $\ast$ denotes the usual convolution product of two scalar-valued distributions$^{[2]}$. Turns out that 
$$\left(\left[\mathscr{D}'(\mathbb{R}^2)\right]^2, +, \star, \cdot_{\mathbb{K}}\right)$$
is an associative and commutative algebra with unity 
$$\begin{bmatrix} 
\delta \\ 
0
\end{bmatrix}.
$$
Moreover, an explicit computation reveals that 
$$
L\left( \begin{bmatrix} 
u_1 \\ v_1\end{bmatrix} \star 
\begin{bmatrix} 
u_2 \\ v_2 \end{bmatrix}
\right) 
=
L\left( \begin{bmatrix} 
u_1 \\ v_1\end{bmatrix}
\right) \star 
\begin{bmatrix} 
u_2 \\ v_2 \end{bmatrix}.
$$
Therefore, if we are able to solve the equation 
$$
L\left( \begin{bmatrix} 
u_E \\ v_E\end{bmatrix} \right) = \begin{bmatrix} \delta \\ 0 \end{bmatrix}, $$
we can later solve the most general inhomogeneous equation
$$
L\left( \begin{bmatrix} 
u \\ v\end{bmatrix} \right) = \begin{bmatrix} f \\ g \end{bmatrix}
$$
by setting $\begin{bmatrix} 
u \\ v\end{bmatrix} = \begin{bmatrix} 
u_E \\ v_E\end{bmatrix}\star \begin{bmatrix} f \\ g \end{bmatrix}$. We can therefore righteously say that $\begin{bmatrix} u_E \\ v_E\end{bmatrix}$ is a fundamental solution to our linear operator $L$. 

$^{[1]}$. The genesis of this formula is the following. You can identify any vector $[x_0, x_1]\in \mathbb{K}^2$ with a $\mathbb{K}$-valued function defined on the cyclic group $\mathbb{Z}/2\mathbb{Z}$. Those groups enjoy a built-in notion of convolution, given by
$$(F\ast G)(n)=\sum_{m\in \mathbb{Z}/2\mathbb{Z}} F(n-m)G(m).$$
Combining this formula with the distributional convolution product you get (1).
$^{[2]}$. Of course we are ignoring here some technicalities regarding the fact that the convolution of two arbitrary distribution is not defined: to be more precise we should restrict our attention to some subspace of distributions having some decay condition at infinity, so that $\ast$ is well-defined and associative.
