@user64494 Not relevant since we do not have absolute convergent series here.

@user64494 Theorem: the sum of two convergent series is convergent. Theorem: if a(k) is alternating and their absolute values converges monotonically to 0, the series with terms a(k) converges. Combine these two.

@user64494 Both series converge according to the Leibniz test for alternating series, so the split into these two series looks fine.

Thanks, and nice! I was confused there for a while!

I fail to see how your formula gives the derivative. In the Cauchy formula for the derivative the $z-z_0$ should be squared, and there should be $2\pi i$ instead of $2\pi$. What do I miss?

This is very useful. Might I ask here if this method also works if the potential is a piecewise polynomial (say $p_1$ for negative $x$ and $p_2$ for positive $x$)? What about if we have an operator on a halfline? Is it possible to make it work? I could ask a new question, if that is the right way to go.

Thank you very much for your answer. I got your code working, and this is exactly in the line I was looking for! I will surely be able to tune it to fit my needs. Also, I'm happy you answered before the bounty expired, so that I can award it to you.

@HenrikSchumacher Thank you for your comments, and many relevant questions. I am not really sure that there is a largest domain. Looking at the situation for different $\Omega_\alpha=\{(x,y)\in D(0,1)~|~y<\alpha\}$ somehow suggests it, but of course, if you have several components, that could be a different story. I would be happy to find a method that somehow starts with the domain in the second example in the question, and somehow pushes the boundary (I know that is not uniquely defined) as far as possible while keeping the negative solution. A bit vague, I know. By largest I mean in area.

@Mathe172 Thank you for your suggestion. I don't see exactly how that applies. Do you mean that I should discretize the disc and then solve the Poisson equation on parts of the mesh in a trial and error way?

I'm sorry, but I don't know how the eigenvalue should behave for larger $h$. It should behave like $$\lambda_1(h)\sim e^{-1/(2h)}(1+O(h))$$ as $h\to 0$ for the disk. At the moment I am clueless.

Thank, you! I just played with that (found it in "MeshRefinementFunction" at the documentation of ToElementMesh. It gets better, but not close to the expected values. I'm afraid it really has to be localized too much for what is practically possible.