You don't need to solve four different boundary value problems: two are enough.
- replace the conditions on vertical sides with homogeneous ones;
- replace the conditions on horizontal sides with homogeneous ones.
Then add. I'll work through part 1.
The reason to have homogeneous BCs on a pair of opposite sides is that we can build them into the solution by using a basis of eigenfunctions of the second derivative operator. In problem 1, we work with eigenfunctions of $\partial^2/\partial x^2$ on $(-a,a)$.
The eigenfunctions for Neumann (or Dirichlet) conditions are trigonometric: sines, cosines, and their linear combinations. In this case, we are to expand $-x/2$, which is an odd function. So we will not need cosines.
Sine function $\sin \beta x$ satisfies the boundary condition $(\sin \beta x)'\bigg|_{x=a}=0$ when $\beta = \frac{2n+1}{2a}\pi$ for some integer $n\ge 0$. I write $\beta_n = \frac{2n+1}{2a}\pi$ below.
Having taken care of the homogeneous boundary conditions, we turn to the PDE itself. The building blocks (separated solutions) for the Laplace equation are of the form "(hyperbolic trig) times (trig)". Namely, any function of the form
$$
u(x,y)=\sum_{n=0}^\infty (A_n\cosh \beta_n y+B_n\sinh\beta_n y) \sin\beta_n x \tag{1}
$$
solves the PDE.
Having taken care of the PDE, we turn to the remaining, inhomogeneous boundary conditions.
The goal is to choose $A_n,B_n$ so that $u$ satisfies them. Observing that $\partial u/\partial y$ is subjected to the same condition on top and bottom sides, we conclude that it should be an even function of $y$, hence $u$ must be an odd function of $y$. This simplifies (1) to
$$
u(x,y)=\sum_{n=0}^\infty B_n\sinh\beta_n y \, \sin\beta_n x
$$
with the derivative
$$
u_y(x,y)=\sum_{n=0}^\infty \beta_n B_n\cosh\beta_n y \, \sin\beta_n x
\tag{2}$$
When $y=b$, the sum (2) must match $-x/2$.
The coefficients of $-x/2$ in the basis $\{\sin \beta_n x\}$ are
$$
\frac{1}{a}\int_{-a}^a -\frac{x}{2}\sin \beta_n x\,dx
= \frac{(-1)^{n+1}}{a \beta_n^2}
$$
unless I messed up the integration. Hence, $\beta_n B_n = (-1)^{n+1}/(a \beta_n^2)$, which yields $B_n= (-1)^{n+1}/(a\beta_n^3)$. Solution:
$$
u(x,y)= \sum_{n=0}^\infty (-1)^{n+1}\frac{8a^2}{(2n+1)^3\pi^3} \sinh \frac{2n+1}{2a}\pi y \, \sin \frac{2n+1}{2a}\pi x
\tag{3}$$
This completes part 1.
Luckily for you, Part 2 is so similar that relabeling $x\leftrightarrow y$, $a\leftrightarrow b$ is all you need.
The solution is
$$
v(x,y)= \sum_{n=0}^\infty (-1)^{n+1}\frac{8b^2}{(2n+1)^3\pi^3} \sinh \frac{2n+1}{2b}\pi x \, \sin \frac{2n+1}{2b}\pi y
\tag{4}$$
The sum $u+v$ is $\phi$ that you wanted.
Some remarks.
The advantage of homogeneous BCs restrict the space of solutions to a linear subspace. When
we look for separated solutions $X_n(x)Y_n(y)$, we can use $X_n$ (or $Y_n$, in part 2) which satisfy the homogeneous BCs; this way, and sum $\sum X_nY_n$ will satisfy them as well. Simply put, $0$ is the only number which we can add to itself, as many times as we want, without changing it.
Summary of the process:
- Two homogeneous BCs are built into the choice of eigenfunction basis.
- The PDE is built into the form of $X_nY_n$ products
- The inhomogeneous BCs are matched by coefficients.
Best Answer
I'm afraid that what I'm about to say is not quite you what to hear: your Neumann PDE does not have a solution for arbitrary choices of $g$.
To see why, let's compute the integral of $g$ over the surface $\Gamma$: $$ \int_\Gamma g \ dS = \int_\Gamma \frac{\partial u}{\partial n} \ dS=\int_\Omega \nabla^2 u \ dV =\int_\Omega 0\ dV = 0.$$ Therefore, the average value of $g$ on the boundary surface $S$ must be zero. If $g$ fails to satisfy this condition, your PDE cannot be solved.
For a similar reason, your definition of the Green's function must be modified. Indeed, if we integrate $\frac{\partial G}{\partial n}$ over the surface $\Gamma$, we find that $$ \int_\Gamma \frac{\partial G(\vec r, \vec r')}{\partial n_{\vec r}}dS(\vec r) = \int_\Omega \nabla_{\vec r}^2 G(\vec r,\vec r') \ dV(\vec r) = \int_V \delta(\vec r-\vec r') \ dV(\vec r) = 1.$$ So it is inconsistent to set $\frac{\partial G}{\partial n}$ equal to zero on the boundary $\Gamma$.
Fortunately, all is not lost! We can redefine the Green's function $G$ so that it satisfies $$ \begin{cases} \nabla_{\vec r}^2 G(\vec r,\vec r') = \delta(\vec r - \vec r') &{\rm \ \ on \ \ } \Omega , \\ \frac{\partial G(\vec r, \vec r')}{\partial n_{\vec r}}=\frac 1 A & {\rm \ \ on \ \ } \Gamma, \end{cases} $$ where $A = \int_\Gamma dS$ is the area of the boundary surface $\Gamma$.
Now, Green's identity states that \begin{multline}\int_\Omega \left( u(\vec r) \nabla_{\vec r}^2 G(\vec r , \vec r') - G(\vec r,\vec r') \nabla_{\vec r} u(\vec r) \right) \ dV(\vec r) \\ = \int_\Gamma \left( u(\vec r) \frac{\partial G(\vec r,\vec r')}{\partial n_{\vec r}} - G(\vec r,\vec r') \frac{\partial u(\vec r)}{\partial n_{\vec r}}\right) dS(\vec r).\end{multline} Plugging in my new definition of $G(\vec r, \vec r')$, we see that any $u(\vec r)$ satisfying your PDE must satisfy we can immediately see that any solution to the PDE must satisfy $$ u(\vec r') = - \int_\Gamma G(\vec r,\vec r')g(\vec r) \ dS(\vec r) + c, \ \ \ \ \ \ (\ast )$$ where the constant $c$ is equal to $\frac 1 A \int_\Gamma u(\vec r) dS(\vec r)$, the average value of $u$ on $\Gamma$.
Having redefined the Green's function, I'll give you an explicit expression in the case where $\Omega$ is a two-dimensional circular disk of radius $1$. Here it is: $$ G(\vec r, \vec r') = \tfrac 1 {2\pi} \ln | \vec r -\vec r' | + \tfrac 1 {2\pi} \ln |\vec r - \vec r''|,$$ where $\vec r'' = \vec r'/|\vec r'|^2$ is the image of $r'$ under an inversion about the unit circle. It is similar to the Dirichlet Green's function, expect that we have a plus sign in front of the "image" term instead of a minus sign.
I believe that if you plug this Green's function into $(\ast )$ (with the constant $c$ chosen arbitrarily), you do indeed get a solution to your PDE, provided that $g$ obeys the consistency condition $\int_\Gamma g \ dS = 0$. This is stated in these lecture notes from my university, and also in Riley, Hobson and Bence, though I would love to see a rigorous proof! I wonder if anyone knows a good reference?