I could not accept my failure and kept improving.
I noticed the problem increases the closer we get to the vertical axis, so i used that instead a damping constant.
It's pretty good now, even in a hard case like this, it nicely converges:

Though, to make the printed error zero, i'd need 42 iterations. It would need more work to map the error better.
So it surely is no great solution, but i found no cases where it does not work.

My improvement form before turned out a case of ‘make it seemingly better by making it actually worse’.
I'll completely replace the code from upper post - too much changes to mention them all.

I could not accept a working but slow solution.

Now it's really good, and i can recommend using it. I think 4 iterations are enough in any case, which was may goal. :D

Instead projecting to ellipse tangent line, i constructed a inner circle at the intersection of ellipse normal and up vector. It's neat. No magic constants, no damping, no thresholds.

I'm happy, and guess it's not much worse than Eberlys maybe, which i still should take a look at now…

Eberlys paper has even code.
And the code looks faster than mine, i'd guess something like 5-10x faster.
But he also states, for floats, his maximum iteration count is 149.

So i did this for a comparison, of which i don't know how fair it is:

					float move = vec2(ellipsePoint - oldEllipsePoint).Length(); // travelled distance withi current iteration
					ImGui::TextColored(ImVec4(R,G,B,1), "move %f", move);
					if (move < std::numeric_limits<float>::epsilon())
						break;				

It shows i need only 6 iterations in difficult cases to reach precision limits.

Thus i may even win over Eberlys method? How cool is that?

Or just Blasphemy? I'd reference him. He was here, but can't find his username.

Thaumaturge said:
If (and as far as) I understand your method correctly, it does seem pretty good, and should produce good results, I imagine!

This is how it works:

You can model an ellipse with circles along it's major axis.
So i get an approximate close point, construct the circle touching that point, calculate exact closest point on the circle, project it to the ellipse and proceed.
The projection errors reduce quickly, so it's fast, and i don't see any failure cases. (Though i have not worked out special cases causing divisions by zero, for example)

Problem with the former method of projecting the query point to the tangent was that this projection was often far away from the current guess on the ellipse, causing oscillations.
With the circles this can't happen because they are not infinite, opposed to a tangent.

Eberlys method surely is more elegant, as he probably constructs a equation to minimize distance, then working on just that ignoring geometry. The iterations are just 5 lines of code of root finding without a need to access memory.
And i guess it works well with much less than 150 iterations, so i really don't know which is faster. Accuracy is probably the same.