I need to solve Laplace equation for a lens (drop on a surface) with constant boundary condition on the top surface and zero on the bottom surface (this is the simplest case, not the only one I considered).
The region and relevant variables and constants are presented in the picture:
Sperical coordinates seem most fitting, so the equation becomes:
$$\frac{1}{r^2}\frac{\partial}{\partial r} \left( r^2 \frac{\partial f}{\partial r} \right)+\frac{1}{r^2 \sin (\theta)}\frac{\partial}{\partial \theta} \left( \sin (\theta) \frac{\partial f}{\partial \theta} \right)=0$$
$$f(r_0 , \theta)=\begin{cases}f_0 & \cos (\theta) > \frac{l}{r_0} \\0 & \cos (\theta) \leq \frac{l}{r_0} \end{cases} $$
$$f\left(r , \arccos \frac{l}{r} \right)=f\left( \frac{l}{\cos (\theta)} , \theta \right)=0 $$
The only case I can work with is the case of a half-sphere, because then $l=0$ and the conditions allow me to build the solution in the form of:
$$f=\sum_{k=0}^{\infty}C_k r^k P_k (\cos (\theta))$$
I solved this equation numerically in Mathematica, but I would like to know how to obtain some sort of approximate solution, like the series above. I researched this problem for months, and was not able to find anything useful.
What I have done so far:
1) I considered an easier problem with the zero boundary condition at the center of the drop, or on the surface of a sphere of a lesser radius. However, these solutions give different size and shape dependences than the numerical solution.
2) I solved this problem in 2D for horizontal cylindrical segment using Mobius transform, and obtained an accurate solution. It is still wrong for the original problem though, and I doubt it's useful for approximation.

The solution: $$f=\sum_{k=0}^{N}C_k r^k P_k (\cos (\theta)) + \sum_{k=0}^{N}D_k r^{-(k+1)} P_k (\cos (\theta)), \;\; N\rightarrow \infty$$ is valid for all axial simmetries.
What you need to do is to find the $C$ and $D$ coefficients that make the solution null on your boundary.
So, for example, as the legendre polynomials are all $1$ for $\cos (\theta)=1$: $$f(l,0)=\sum_{k=0}^{N}C_k l^k + \sum_{k=0}^{N}D_k l^{-(k+1)}=0$$
Using enough of this kind of expressions you can find always better approximations to the real solutions: choose an $N$ large enough, evaluate the function in $2N$ points and find the approximate values for $C_k$ and $D_k$.
Take care when angles are near $\cos (\theta)=l/r_0$ as you have a discontinuity on the boundary and this will make the system ill-conditioned. Better not to sample the boundary near those points.