Problem solving PDE with boundary conditions

Question

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,...$.

Answer 1

Something like this should get you started:

(*analyticalSolution[x_,y_]:=...*)
bcs = DirichletCondition[\[Phi][x, y] == analyticalSolution[x, y], 
   True];

sol = NDSolveValue[{Inactive[
      Div][{{Gx, 0}, {0, Gy}}.Inactive[Grad][\[Phi][x, y], {x, 
         y}] + {{-Gx y, Gy x}}, {x, y}] == 0(*,bcs*)}, \[Phi], {x, -a,
    a}, {y, -b, b}]