# Problem solving PDE with boundary conditions

I have the following PDE: $$\frac{\partial }{\partial x}\left(G_x \left(\frac{\partial \phi (x,y)}{\partial x}-y\right)\right)+\frac{\partial }{\partial y}\left(G_y \left(\frac{\partial \phi (x,y)}{\partial y}+x\right)\right)=0$$ with BCs $$\frac{\partial \phi (x,y)}{\partial x}-y=0$$ at $$x=\pm a$$ and $$\frac{\partial \phi (x,y)}{\partial y}+x=0$$ at $$y=\pm b$$, $$G_x$$ and $$G_y$$ are constants.

After a lot of work and help from this community I managed to get an analytical solution. Now I want to confirm this solution using NDSolve so I typed the following code:

(*Main equation*)
eqn[x_,y_]=D[(Gx(D[\[Phi][x,y],x]-y)),x]+D[(Gy(D[\[Phi][x,y],y]+x)),y];
(*BCs*)
bcx[x_,y_]=D[\[Phi][x,y],x]-y;
bcy[x_,y_]=D[\[Phi][x,y],y]+x;
bcs={bcx[-a,y]==0,bcx[a,y]==0,bcy[x,-b]==0,bcy[x,b]==0};
(*Values for numerical solution*)
Gy=41018756.0;Gx=72463203.0;a=0.0025;b=0.0025;

NDSolve[{eqn[x,y]==0,bcs},\[Phi],{x,-0.0025,0.0025},{y,-0.0025,0.0025}]


but I get the message

NDSolve::fembdnl: The dependent variable in -y+(\[Phi]^(1,0))[-0.0025,y]==0 in the boundary condition DirichletCondition[-y+(\[Phi]^(1,0))[-0.0025,y]==0,x==-0.0025] needs to be linear.


I thought that all equations can have a numerical solution but are not guaranteed an analytical one so I'm sure I'm making a mistake but I can't see where.

Side note: here is the analytical solution $$\phi (x,y)=x y-\frac{32 \sqrt{G_y} (-1)^n \sin \left(\frac{1}{2} \pi (2 n+1) x\right) \text{sech}\left(\frac{\pi b \sqrt{G_x} (2 n+1)}{2 \sqrt{G_x}}\right) \sinh \left(\frac{\pi \sqrt{G_x} (2 n+1) y}{2 \sqrt{G_y}}\right)}{\pi ^3 \sqrt{G_x} (2 n+1)^3}$$ where $$n=0,1,2,3,...$$.

(*analyticalSolution[x_,y_]:=...*)