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

We present a numerical approach which allows the solving of Bethe equations whose solutions define the eigenstates of Gaudin models. By focusing on a new 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 Richardsons fermionic pairing model, the central spin model and generalized Dicke model.