The part of the iterative method is not easy to understand.
I found the document on which you rely on the following link:
The latitude of a point on the central meridian in a conformal Gauss projection is found considering its northing as the meridian arc length from the point to the equator, multiplied by the scale (k_0) of the projection.
The calculation of the meridian arc length (S) is not being performed by an iterative method but rather as the development in a truncated series of the ellipse integral.
What is not clear is whether we want to calculate the latitude of the point on the central meridian (perpendicular projection of the point in question towards the central meridian) using the northing of the point in question.
Without a doubt, if that were done and the latitude of the point was calculated from there, then its northing was calculated, there would be a difference with the original northing of the point. It is not understood if it is there where it is intended to iterate.
Iterative methods are often used to find ellipsoidal distances between points through spherical approximations, since the parameters of the ellipse containing the points are not known. The solution of the ellipse integral, when its parameters are known, is usually found as a series expansion, truncating the series (usually in the fourth term) rather than by numerical methods.
About the implementation in the library to which you refer, it seems to follow the formulas provided in the wikipedia article:
All the links referred to in that article are of great interest, beginning with the original development of Krüger (1912):
An interesting analysis and further development is what Karney (2011) performs:
And the simplified formulas provided in the article appear to be those developed by Kawase (2012):