$u_{xx}+u_{yy}=0$ on $0
$u_{xx}+u_{yy}=0$ on $0<x<a,0<y<b$

$u_{x}(0,y)=-a$, $u_x(a,y)=0$


$u_y(x,0)=b$, $u_y(x,b)=0$

I separated this into $u=v+w$ where

$v_{xx}+v_{yy}=0$ on $0<x<a,0<y<b$


$v_{x}(0,y)=-a$, $v_x(a,y)=0$


$v_y(x,0)=0$, $v_y(x,b)=0$

and

$w_{xx}+w_{yy}=0$ on $0<x<a,0<y<b$


$w_{x}(0,y)=0$, $w_x(a,y)=0$


$w_y(x,0)=b$, $w_y(x,b)=0$

Solving for $v$
I used separation of variables to get $v=X(x)Y(y)$
And $\frac{X''}{X}+\frac{Y''}{Y}=-\lambda$ where $\lambda\geq 0$
Which gives me ODEs, $X''+\lambda X=0$, $Y''+\lambda Y=0$
$Y(y)=A\cos(\sqrt{\lambda} y)+B\sin(\sqrt{\lambda}y)$
Initial condition $v_y(x,0)=0$ gives me $B=0$
So I get $Y=A\cos(\sqrt{\lambda}y)$
And then the condition $Y'(b)=0$ gives me that $\sqrt{\lambda}=\frac{n\pi}{b}$
So $Y=A\cos(\frac{n\pi}{b})$
And solving for $X$ I get $X=De^{\frac{n\pi}{b}x}+Ce^{-\frac{n\pi}{b}x}$
I end up with this general solution for $v(x,y)= \sum_{n=1}^\infty \cos(\frac{n\pi}{b}y)[A_ne^{\frac{n\pi}{b}x}+B_ne^{-\frac{n\pi}{b}x}]$
My issue now is when I try to use $v_x(0,y)=-a$ and $v_y(a,y)=0$ to solve for $A_n, B_n$
I get $v_x(0,y)=-a=\sum_{n=1}^\infty \cos(\frac{n\pi}{b}y) [A_n \frac{n\pi}{b}-\frac{n\pi}{b}B_n]$
and $v(a,y)=0= \sum_{n=1}^\infty \cos(\frac{n\pi}{b}y)[A_n\frac{n\pi}{b}e^{\frac{n\pi}{b}a}-B_n\frac{n\pi}{b}e^{-\frac{n\pi}{b}a}]$
Which do not seem to help me solve for $A_n$ and $B_n$.
 A: Note that your general solution for the Neumann problem should also include the $n=0$ (constant term).  Then in principle, $A_n$ and $B_n$ would come from matching the coefficients in your solution to the Fourier cosine coefficients of the boundary conditions.
The more serious challenge is that your decomposition $u = v+w$ will not work, due to the solvability condition for the Neumann problem for Laplace's equation; the total flux $\int_{\partial\Omega}\frac{\partial u}{\partial n} dl$ into the domain must be zero.  While the boundary conditions for $u$ satisfy this condition ($ab - ba = 0$), the ones for $v$ and $w$ do not.
There is more than one way around this.  One could be to do the same decomposition $u = v+w$, with your boundary conditions, but with $u$ and $v$ instead satisfying
$$
v_{xx} + v_{yy} = C, \qquad w_{xx} + w_{yy} = -C
$$
where $C$ is a constant (only a unique value of $C$ will allow for a solution to each problem).  Each subproblem will be greatly simplified by the observation that the solution for $v$ is independent of $y$, and the solution for $w$ is independent of $x$.  In the end you will get a quite simple polynomial solution in $x$ and $y$.
