In this work, we generalize the numerical approach to Gaudin models developed earlier by us to degenerate systems showing that their treatment is surprisingly convenient from a numerical point of view. In fact, high degeneracies not only reduce the number of relevant states in the Hilbert space by a non negligible fraction, they also allow to write the relevant equations in the form of sparse matrix equations. Moreover, we introduce a new inversion method based on a basis of barycentric polynomials which leads to a more stable and efficient root extraction which most importantly avoids the necessity of working with arbitrary precision. As an example we show the results of our procedure applied to the Richardson model on a square lattice.