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
Separate variables as $u(x_1,x_2) = F(x_1)G(x_2)$. Laplace's equation becomes (upon division by $FG$): \begin{equation} \frac{F''}{F} =- \frac{G''}{G} = \lambda \in \mathbb{R}. \end{equation} This gives two possibilities depending on the sign of $\lambda$: oscillatory or exponential behaviour. From the data, $F(x_1)G'(0) = \frac{1}{n} \sin(nx_1)$, so we may absorb constants and WLOG take $G'(0)= 1$, $F(x_1) = \frac{1}{n} \sin(nx_1) $.
This then determines the value $\lambda = -1/n^2$. Then, you obtain a general solution for $G$: \begin{equation} G(x_2) = A \cosh(nx_2) + B \sinh(nx_2). \end{equation} You have two boundary conditions $G(0) = 0$, $G'(0) = 1$ which tell you that $A=0$, $B = 1/n$ as needed.
As $n \to \infty$ the solution blows up ( take $x_2 \neq 0, x_1 = \pi / 2n$). In general, Cauchy problems are only posed for hyperbolic equations and Laplace's equation is an elliptic equation. In this case, as $n \to \infty$ the initial data approaches $u=0, u_{x_2}=0$ at $x_2 = 0$, which gives rise to the solution $\tilde{u} \equiv 0$. However, $|u-\tilde u|$ can be made arbitrarily large as $n \to \infty$ (take $x_1 = \pi /2 n$). Thus, there is no continuous dependence on the data (perturbing one boundary condition gives a very different solution), so the problem is not well-posed.