  The study of the electrodynamics of static, axisymmetric and force-free Kerr
magnetospheres relies vastly on solutions of the so called relativistic
Grad-Shafranov equation (GSE). Different numerical approaches to the solution
of the GSE have been introduced in the literature, but none of them has been
fully assessed from the numerical point of view in terms of efficiency and
quality of the solutions found. We present a generalization of these algorithms
and give detailed background on the algorithmic implementation. We assess the
numerical stability of the implemented algorithms and quantify the convergence
of the presented methodology for the most established setups (split-monopole,
paraboloidal, BH-disk, uniform).
