Method for solving collection of simple PDEs How would you go about evaluating the following collection of simple PDEs: 
$$\frac{\partial A_3}{\partial y} - \frac{\partial A_2}{\partial z} = yz$$
$$\frac{\partial A_1}{\partial z} - \frac{\partial A_3}{\partial x} = xz $$
$$\frac{\partial A_2}{\partial x} - \frac{\partial A_3}{\partial y} =xy $$
I'm trying to solve for $A_1(x,y,z), A_2(x,y,z), A_3(x,y,z)$ in a methodical way. If you have a good reference to notes that would be great. Thanks for any assistance. 
 A: Disclaimer 
I have no idea what I am doing. I would love to investigate this direction more deeply and to look at some actual theory. The following is basically me trying to go ahead with the calculations; as such I am sure there will be lots of errors / unstated assumptions.
The Linear Algebra Point of View
Let's view this as a linear algebra problem.We isolate the partial derivatives, that we now see as a linear operators. We get 
$$Ba = c$$
where $$a = \begin{pmatrix} A_1 \\ A_2 \\ A_3 \end{pmatrix} \ \ \ \ \ \  c = \begin{pmatrix} yz\\ xz \\ xy \end{pmatrix} \ \ \ \ \text{       and } B = \begin{pmatrix} 0 & -\partial_z & \partial y  \\ \partial_z &0 &-\partial_x  \\ 0 & \partial_x & -\partial_y \end{pmatrix}$$
Where $B$ represents a linear transformation on the space of (differentiable) functions $f(x,y,z)$. To make life easier we can take $f$ to be smooth, or at least $C^2$ (basically whatever order we need to make the following calculations work). In this way we are possibly losing some (non-smooth) solutions but let's forget about this for now. 
The next logical step is that we want to invert the matrix $B$ and multiply it by $c$. What do we need to be able to invert a matrix? We need the notion of multiplication, addition, and inverse of an element ("division"). Addition and multiplication are easy: in fact we can see $$(\partial_z + \partial_y) f = \partial_z f + \partial_y f$$
$$(\partial_z\partial_y)f = \partial_z(\partial_y f) = \partial_y(\partial_z f)$$
Inverse of an element (a derivative) must be "the integral" of some sort; more generally it's a solution to a partial differential equation. So for example we write 
$$\partial_z^{-1} f = \int f dz$$ 
$$(\partial_x - \partial_z)^{-1} f = g$$ where $g$ is "one" solution to the PDE
$$\partial_x g - \partial_z g = f$$
Notice how $g$ can be interpreted as a whole class of functions, as there is no single solution to that PDE without boundary conditions. But let's keep going for now and see where we end up.
For notation sake we will also write $\frac 1{\partial z}$ instead of $\partial_z^{-1}$ when it's convenient to do so.
The inverse of $B$ is then given by 
$$\large B^{-1} = \begin{pmatrix} \frac{\partial_x \partial_x}{ \partial_y \partial_z(\partial_x - \partial_z)} & \frac{1}{\partial_z} & \frac{\partial_x}{\partial_y(\partial_x - \partial_z)}  \\ \frac{1}{\partial_x - \partial_z} & 0 & \frac{1}{\partial_x - \partial_z} \\ \frac{\partial_x}{\partial_y(\partial_x - \partial_z)} & 0 & \frac{\partial_z}{\partial_y(\partial_x - \partial_z)} \end{pmatrix}$$
Therefore the solution is given by 
$$\begin{pmatrix} A_1 \\ A_2 \\ A_3 \end{pmatrix} = B^{-1} \begin{pmatrix} yz  \\ xz  \\ xy \end{pmatrix}$$
Notice that we still have to solve PDEs, but now each PDE is not "connected" to the other ones. It's not a system of PDEs anymore but just multiple PDEs to solve.
Carrying on with the computation, we find that $A_1$ is given by 
$$A_1 = \frac{\partial_x \partial_x}{ \partial_y \partial_z(\partial_x - \partial_z)} (yz) +  \frac{1}{\partial_z}  (xz)  +  \frac{\partial_x}{\partial_y(\partial_x - \partial_z)}(xy)$$
Now we need to prove the following property about the "inverse" symbol we defined above, namely that $$\frac 1{\partial_y(\partial_x - \partial_z)} = \frac 1{\partial_y}\cdot \frac 1{\partial_x - \partial_z}$$
This is very intuitive and means that we can solve solve the "bigger" pde solving a smaller one and putting that general solution as "known term" for the other one.  Proof in Appendix 1
Using this property we get that $$A_1 = \frac 0{\partial_y \partial_z (\partial_x - \partial_z)} + xz^2/2 + \frac{y^2/2}{\partial_x - \partial_z}$$
For the first term, we know that the general solution to $\partial_y \partial_z f = 0$ is $f = L(y) + H(z)$ for some functions $L, H$. Hence we can write 
$$A_1 = \frac{L(y) + H(z) + y^2/2}{\partial_x - \partial_z} + xz^2/2$$.
Doing the same for the other two we get 
$$A_2 = \frac{yz + xy}{\partial_x - \partial_z}$$
$$A_3 = \frac{D(x,z)}{\partial_x - \partial_z}$$
for some function $D(x,z)$ that does not depend on $y$. 
But the functions $L,H,D$ cannot be chosen independently [1]. In fact substituting $A_1, A_2, A_3$ in the original PDEs we were trying to solve we get that 
$$\frac {-y}{\partial_x - \partial_z} = yz \iff -y = (\partial_x - \partial_z)(yz)$$ This is always true, so no conditions from this equation.
$$\frac{H'(z)}{\partial_x - \partial_z} + xz - \frac{\partial_x D}{\partial_x - \partial_z} = xz \iff H'(z) = \partial_x D(x,z)$$
Finally the third:
$$\frac{y}{\partial_x - \partial_z} = xy$$ which is true, so no conditions. All in all we finally have 

$$A_1 = \frac{L(y) + H(z) + y^2/2}{\partial_x - \partial_z} + xz^2/2$$
  $$A_2 = \frac{yz + xy}{\partial_x - \partial_z}$$ $$A_3 =
 \frac{D(x,z)}{\partial_x - \partial_z}$$
with $H'(z) = \partial_x D(x,z)$

Note how these are not closed form solution because we still have to solve the PDEs where the partial derivatives are at the denominator.
Now choosing for example $D = H = L = 0$, we get 
$$A_1 = \frac x2(z^2 + y^2)$$
$$A_2 = \frac y2(x^2 - z^2)$$
$$A_3 = 0$$
(I skipped over another problem; namely how to choose the solution from the last PDEs (even when setting $L, H, D$ to $0$). The solution I chose works but there are other choices that don't, and I am not sure why. Probably there are more conditions to put but not sure where they would come from. )
[1] For some reason.. I thought that they could. Apparently not, I'll have to think more over this.
Appendix 1
We have  $$\frac 1{\partial_y(\partial_x - \partial_z)}f = g$$ where $g$ is the solution to $$\partial_y(\partial_x - \partial_z) g = f$$. Calling $(\partial_x - \partial_z) g  = h$ we then have $\partial_y h = f \implies h = \int f dy \implies (\partial_x - \partial_z)g = \int f dy$
which implies that 
$$g =  \frac 1{\partial_x - \partial_z}\left(\frac 1{\partial_y}f\right)$$
