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 Richardson's fermionic pairing model, the central spin model and generalized Dicke model.The correspondence between Bethe Ansatz / Ordinary Differential Equation (ODE) has been found some time ago 1 on the basis of similarity between T-Q system of the Bethe Ansatz and functional relations including spectral determinants on the ODE side. For extensive review we refer to 2 . The correspondence is related to the Langlands correspondence where the ODE part has to do with socalled Miura opers 3 . It has many interesting links with conformal field theory, quasi-exact solvability, supersymmetry and PT-symmetric quantum mechanics. It was also used it in the context of physics of cold atoms and quantum impurity problem 4 . Here we develop further this remarkable correspondence to facilitate the solving of otherwise difficult-to-solve Gaudin systems and implement it to several physically interesting situations.Indeed, numerous integrable models derived from a generalized Gaudin algebra have been used to describe the properties of physical systems. The Richardson fermionic pairing Hamiltonian 5 has explained most properties of superconducting nanograins 6,7 , numerous studies of the decoherence of a single electron spin trapped in a quantum dot rely on the central spin model 8 while the inhomogeneous Dicke models 9 have been used to study light-matter interaction in many cavity-based systems.The quantum integrability of such systems brings major simplifications to the structure of their eigenstates. In a system containing N degrees of freedom, they can be fully described by a number M of complex parameters we call rapidities. M is here a number of the same order as N despite the exponentially large Hilbert space. The possible values of this set of parameters can be obtained by finding solutions to a set of M coupled non-linear algebraic equation: the Bethe equations. However, analytically finding solutions to the Bethe equations remains impossible except in very specific cases. The problem is still numerically approachable but solving non-linear equations remains a challenge that only iterative methods can tackle.While previous efforts in designing algorithms have shown promising results [10][11][12][13][14], methods which are sufficiently fast and stable for systematically finding a large number of eigenstates remained elusive. As was shown by one of the authors, computing efficiently only a small fraction of the complete Hilbert space can be sufficient to access static 15 , dynamical 16 and even non-equilibrium dyna...