Problem:
Consider the BVP:
\begin{equation}
\label{eq1}
\begin{cases}
\frac{\partial^2\Psi}{\partial\rho^2} + \frac{1}{\rho}\left(\frac{1-h^2\rho^2}{1+h^2\rho^2}\right)\frac{\partial\Psi}{\partial\rho}+\frac{1+h^2\rho^2}{\rho^2}\frac{\partial^2\Psi}{\partial\zeta^2} = \frac{-2hF_0}{1+h^2\rho^2}\,, & \text{on } \Omega\,,\\
\Psi(\rho,\zeta) = 0\,, & \text{on } \partial\Omega\,.
\end{cases}
\end{equation}
where $\Omega$ is the ellipse $\{(x,y)\in\mathbb{R}^2\mid (x-R_0)^2/a^2 + y^2/b^2\leq 1\}$ with $R_0-a>0$ and $\rho,\zeta$ are polar coordinates. I am trying to determine the number of critical points the solution to the above equation has, and the nature of these critical points.
My Attempts:
Numerical Solutions using Mathematica:
The numerical solution to the PDE using NDSolveValue has a single critical point on the line $\zeta= 0$. Moreover, along this line we find that $\frac{\partial\Psi}{\partial\rho}$ seems to be monotonically decreasing. Another point to note is that as $b/a$ increases past $1$, the location of the critical point shifts to the right of $R_0$. Shown below is the solution for $a=1,b=3$.
Analysis
The above question is equivalent to finding the zeroes of $Y =(\frac{\partial\Psi}{\partial\rho},\frac{\partial\Psi}{\partial\zeta})= (Y_1,Y_2)$ satisfying \begin{equation} \label{eq2} \left[\left(\frac{\partial}{\partial\rho},\frac{1+h^2\rho^2}{\rho^2}\frac{\partial}{\partial\zeta}\right) + \frac{1}{\rho}\left(\frac{1-h^2\rho^2}{1+h^2\rho^2}\right)\left(1,0\right)\right]\cdot\left(Y_1,Y_2\right) = \frac{-2hF_0}{1+h^2\rho^2} \end{equation}
By the symmetry of the PDE, we have that $\Psi(\rho, -\zeta) = \Psi(\rho,\zeta)$. Thus $Y_2(\rho, -\zeta) = -Y_2(\rho, \zeta)$ so $Y_2(\rho,0) = 0$. The boundary condition $\Psi = 0$ on $\partial\Omega$ implies that $\nabla\Psi$ points inwards on the boundary, so $Y_1(R_0-a,0)$ and $Y_1(R_0+a,0)$ are both non-zero and have opposite sign. Therefore at some point on the line $\zeta=0$ we have $Y_1(\rho,\zeta)=0$ and so there is at least one critical point on the line $\zeta = 0$.
Looking for solutions along $\zeta=0$
points along $\zeta = 0$:} Along $\zeta=0$, we can consider $Y_1$ as a function of only $\rho$. So the above equation becomes $$ Y_1'(\rho) + \frac{1}{\rho}\left(\frac{1-h^2\rho^2}{1+h^2\rho^2}\right)Y_1(\rho) = \frac{-2hF_0}{1+h^2\rho^2} - \frac{1+h^2\rho^2}{\rho^2}\frac{\partial Y_2}{\partial\zeta}(\rho,0) $$ which implies, $$ (IY_1)' = \frac{-2h\rho F_0}{(1+h^2\rho^2)^2} - \frac{1}{\rho}\frac{\partial Y_2}{\partial\zeta}(\rho,0) $$ where $I = \exp\left(\int_1^\rho\frac{1}{\rho}\left(\frac{1-h^2\rho^2}{1+h^2\rho^2}\right)\,ds\right) = \frac{\rho}{1+h^2\rho^2}$. Since our numerical results suggest that $Y_1$ is monotonically decreasing, one approach would be to show that $(IY_1)' < 0$ everywhere on $\zeta=0$. Unfortunately, $\frac{\partial Y_2}{\partial\zeta}(\rho,0)$ depends on boundary information not just on the line $\zeta = 0$.
Does anyone have any advice for how I can further approach this problem?