Steady free surface (FS) flows can be solved numerically with capturing or fitting methods, the latter being the subject of this paper. Most fitting methods are (pseudo-)transient and thus quite slow for steady flows; the so-called steady iterative method is much faster, but requires a dedicated solver because of the complex FS boundary conditions. The goal is to develop a (currently 2D) fitting method which is fast and can be used with a black box flow solver. Results from a perturbation analysis are used in combination with the IQN-ILS algorithm to construct such a method, applicable to supercritical flows. To tackle this method's scaling problem when the mesh is refined, an extension is proposed which uses a multigrid technique for the s...