In a halfspace, the construction of Neumann-Green function proceeds by reflection similar to the construction of Dirichlet-Green function; the difference is that one uses even reflection instead of odd reflection, thus achieving zero derivative instead of zero value. So, in $n\ge 3$ dimensions the function is (according to one convention)
$$G(x;y) = -\frac{c_n}{|x-y|^{n-2}}-\frac{c_n}{|x-\xi|^{n-2}}$$
where $\xi$ is the reflection of $y$ in the boundary of the halfspace.
On a bounded domain such as a ball, the Neumann-Green function has to be defined a bit differently from the Dirichlet-Green function, because having $\Delta_x G(x;y)=\delta_y$ is incompatible with the Neumann boundary condition. Instead, the requirement on the Laplacian is $$\Delta_x G(x;y)=\delta_y-k,\quad \text{ where } \ k=|\Omega|^{-1}\tag{1}$$ Some authors prefer to put Laplacian on the second variable here; but when normalized by condition $\int G(x,y)\,dx=0$, the function is symmetric.
The requirement (1) makes the integral $\int_\Omega \Delta_x G(x;y)\,dx$ equal to zero, which is consistent with the boundary condition $\frac{\partial }{\partial n_x}G(x;y)=0$ on $\partial\Omega$. Note that the usage of Neumann-Green function is not affected by this $-k$; given any reasonable function $f$ with $\int_\Omega f=0$, we can solve the Neumann problem for the Poisson equation $\Delta u=f$ as
$$
u(x)=\int_\Omega G(x;y) f(y)\,dy
$$
since $$\Delta u(x)=\int_\Omega \Delta_x G(x;y) f(y)\,dy = f(x)- \int_\Omega k f(y)\,dy = f(x)$$
Explicit formulas
The singular part of the Laplacian in (1) is contributed by the fundamental solution; the constant part can be contributed by a quadratic polynomial. The hard part is to find a harmonic function to cancel out the normal derivative. A step toward finding it is to perform even reflection as for half-plane; in two dimensions this is enough but in three dimensions one needs further correction.
In the book Partial Differential Equations by Emmanuele DiBenedetto, pages 106-107, one can find the following formulas for $G(x;y)$ in the ball of radius $R$. (NB: DiBenedetto uses the convention $\Delta_y G(x;y)=-\delta_x+k$.)
Two dimensions
$$G(x;y) = -\frac{1}{2\pi}\left( \ln|\xi-y|\frac{|x|}{R} +\ln|x-y|\right) - \frac{1}{4\pi R^2} |y|^2$$
where $\xi=\frac{R^2}{|x|^2}x$ is the reflection of $x$ across the boundary of the ball.
Three dimensions
$$G(x;y) = \frac{1}{4\pi}\left( \frac{1}{|x-y|} + \frac{R}{|x|} \frac{1}{|\xi-y|}\right)
+ \frac{1}{4\pi}\ln\left( (\xi-y)\cdot \frac{x}{|x|} +|\xi-y|\right)
-\frac{1}{8\pi R^3} |y|^2 $$
Unlike the Dirichlet case, the author does not give a formula for $n$ dimensions, which suggests there isn't a simple one.
I'm not sure in what generality and rigor you are expecting an answer; you may want to look at the theory of pseudodifferential operators on vector bundles. But perhaps the following will be more relevant to you:
It helps to think of $ P $ as an operator with matrix coefficients or equally as a matrix of differential operators.
The identity operator $ \mathcal{I} $ defined by $ \mathcal{I}[v] = v $ has the kernel $ I \; \delta(r) $ where $ I: \mathbb{R}^3 \rightarrow \mathbb{R}^3 $ is the identity matrix,
$$
\mathcal{I}[v](r) = \int \delta(r-r') \; I \cdot v(r') \; d^3r' = \int \delta(r-r') \; I \cdot v(r') \; d^3r' = \int \delta(r-r') \; v(r') \; d^3r' = v(r).
$$
Since $ I $ is a constant, it just factors out of the fourier transform, so that transform of the kernel of $\mathcal{I}$ is just $ I \; e^{i k \cdot r } $.
Let's look at those examples you had. The vector laplacian is
$$
\begin{pmatrix} \Delta & 0 & 0 \\ 0 & \Delta & 0 \\ 0 & 0 & \Delta \end{pmatrix},
$$ and in the fourier basis $$ ( k_x^2 + k_y^2 + k_z^2 ) \; I, $$ so it is easy to invert, and after inverse fourier transform will give you the usual Green's function.
Curl is
$$
\begin{pmatrix} 0 & -\partial_z & \partial_y \\ \partial_z & 0 & - \partial_x \\ -\partial_y & \partial_x & 0\end{pmatrix}.
$$
In fourier basis
$$
C = i \begin{pmatrix} 0 & -k_z & k_y \\ k_z & 0 & -k_x \\ - k_y & k_x & 0 \end{pmatrix}
$$
This one is tricky to invert because it has a kernel (gradients) and cokernel (div of curl is zero). You can check though that
$$
C \cdot \left( \frac{1}{k_x^2 + k_y^2 + k_z^2} C \right)
$$
is the identity when restricted to subspace $ \left\{ v \; \big| \; v \cdot k = 0 \right\} $, which is precisely the subspace of divergence free fields. Inverse fourier transform to obtain
$$
\frac{1}{(x^2+ y^2 + z^2)^{3/2}} \begin{pmatrix} 0 & -z & y \\ z & 0 & -x \\ -y & x & 0 \end{pmatrix},
$$
which is a weird way of writing the kernel
$$
\frac{(r - r')}{|r-r'|^3} \times \cdot
$$
that you will recognize from the Biot-Savart law.
(Sorry for sloppy notation, missing $ -1 $s and $ \pi $s, etc.)
Best Answer
A Greens function $G(x,y)$ is a continuous function such that for each $x$, $$ -\partial_y^2 G(x,y) = \delta_x,\\ G(x,y)|_{y\in\{0,L\}}=0$$
These two conditions are enough to show that $f(x)=\int_0^L g(y)G(x,y) dy$ solves $-f''=g$ with Dirichlet boundary data. If you do not know what the first equation means, you can take it to mean that $\partial_yG$ is a function such that \begin{align} y\neq x &\implies \partial^2_y G(x,y) = 0, \tag{A}\\ \epsilon>0 &\implies \partial_y G(x,x+\epsilon) - \partial_y G(x,x-\epsilon) = -1 \tag{B} \end{align}
(A) means that $\partial_y G$ is locally constant in $y$, away from $y=x$. (B) means that for each $x$, there is a single jump discontinuity in $\partial_y G$ at $y=x$. Hence, for each $x$, $G(x,y)$ is two line segments, one from $y=0$ to $y=x$ and another from $y=x$ to $y=L$. Therefore, there are $A(x),B(x),C(x),D(x)$ such that $$G(x,y) = \begin{cases}A(x)y+B(x) & y<x \\ C(x)y + D(x) & y \ge x \end{cases}$$ with the constraints \begin{align} G(x,0) = 0 &\implies B(x) = 0 \tag{1} \\ G(x,L) = 0 &\implies D(x) = -C(x)L \tag{2} \\ G(x,y) \text{ is $C^0$ at $y=x$} &\implies A(x)x = C(x)(x-L) \tag{3} \\ \partial_y G(x,y+) - \partial_y G(x,y-) = -1 &\implies C(x) = A(x)-1 \tag{4} \end{align}
(4) into (3) gives $$A(x)x=(A(x)-1)(x-L)=A(x)x-x+L-LA(x) \\ \implies A(x)=-(x-L)/L $$ in summary \begin{align} A(x)&=-(x-L)/L \\ B(x) &= 0 \\ C(x) &=-x/L \\ D(x) &= x \end{align} giving the formula.