Specifying boundary conditions for unusual heat-like equation

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}]
• I don't want to solve them as a series of ODE's for 2 reasons: (a) I need to perform a variety of numerical integration and differentiation on F with wrt to y once I get the solution, (b) I will eventually add terms to the model that will have y-derivatives and integrals of F over y. – kotozna Nov 7 '17 at 22:02
• Anyways, I would classify −B[y] K[x] as load, not as boundary condition. – Henrik Schumacher Nov 7 '17 at 22:24
• I tried to simplify my equation when I presented it here, seems I made it too simple. Edited question now contains extra term in y. – kotozna Nov 7 '17 at 22:27
• The sign of $A[x]$ is crucial: It decides if you have an elliptic or hyperbolic equation. Are you sure that you specified it correctly? In case of $A[x]>0$, this should be a hyperbolic equation (similar to the wave equation). In the case $A[x]<0$, it is elliptic. So what should it be? – Henrik Schumacher Nov 7 '17 at 22:38
• @HenrikSchumacher In my real case the sign of $A[x]$ can change. For now what I specified above is fine, I will ask a new question if I have trouble for the other cases. – kotozna Nov 8 '17 at 1:03