from mpmath import mp
mp.dps = 50
k = 3
psi = lambda x: x*(1-mp.e**(-x))/(1-mp.e**(-x)-x*mp.e**(-x))
xi = mp.findroot(lambda x: psi(x) - k, k)
print(xi)Each entry is computed in ball arithmetic with 64 guard bits beyond the 100 digits written. The root $\xi_k$ is enclosed by bisection on the sign of $\psi(x)-k$, down to a bracket of half-width $10^{-106}$, with the signs at the two ends checked in ball arithmetic; that bracket is the stored ball.
The generator compares every row with an mpmath solve at 120 digits, and compares the derived thresholds for $3\leq k\leq7$ with the ten-decimal table in [2].