Specifying boundary conditions for unusual heat-like equation

Question

I am trying to numerically solve the following equation which arises in problems associated with heat flow:

$y^{-3} F[x,y]-\partial_x (y^3 \partial_xF[x,y])+y \partial_x(A[x]y^3 \partial_yF[x,y])=-B[y]K[x]$

where the functions $B[y]$, $A[x]$ and $K[x]$ are known, for example I might use:

$B[y]=y^2 \exp(-y^2)$

and

$K[x]=G[x]^{3/2} \partial_xG[x]$

end

$A[x]=-\partial_xG[x]$

where $G[x]=1+(1+\operatorname{Tanh}[4(x-1/2)])$.

This is a 2D PDE. The values of the function F are not known at the boundaries $x=0$ and $x=L$, but their derivatives are:

$\partial_x F[x,y]=0$ at $x=0$

$\partial_x F[x,y]=0$ at $x=L$

The other BC's are Dirichlet:

$F[x,y=0]=0$

Strictly the other boundary condition for F is only known ay y=infinity. However, I expect the function F to decay strongly with y to zero at infinity, and will set it approximately equal to zero at some relatively large value of y (=Y):

$F[x,y=Y]=0$

First question is what type of solution method can I use for this equation? i.e. is "NDSolve" appropriate? If so how do I implement these boundary conditions. I am confused as to whether the known functions $B[y]$ and $K[y]$ need to be treated "as part of the equation" or as part of the boundary conditions instead. I question using NDSolve this because the boundaries in $x$ are of the Neumann type, and according to the Mathematica documentation, if you have a problem involving Neumann BC's, then you should use finite elements, which I believe requires you to use "NDSolveValue" - is that correct?

Mathematica code for my attempt:

G = 1 + (1 + Tanh[4*(x - 0.5)])

A = (-G^(3/2))*D[G, x]

eqn := F[x, y]/y^3 - D[y^3*D[F[x, y], x], x] + 
   y*D[A*y^3*D[F[x, y], y], x] == (-y^2)*Exp[-y^2]*D[G, x]

sol3 = NDSolve[eqn, F, {x, 0, 1}, {y, 0, 6}]