Gaudin models solver based on the correspondence between Bethe ansatz and ordinary differential equations

We present a numerical approach which allows the solving of Bethe equations whose solutions define the eigenstates of Gaudin models. By focusing on a different set of variables, the canceling divergences which occur for certain values of the coupling strength no longer appear explicitly. The problem is thus reduced to a set of quadratic algebraic equations. The required inverse transformation can then be realized using only linear operations and a standard polynomial root-finding algorithm. The method is applied to Richardson's fermionic pairing model, the central spin model, and the generalized Dicke model.