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/(k*(1-mp.e**(-xi))**(k-1)))
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$, where $\psi(x)=x(1-e^{-x})/(1-e^{-x}-xe^{-x})$, down to a bracket of half-width $10^{-106}$, with the signs at the two ends checked in ball arithmetic; the value of $x/(k(1-e^{-x})^{k-1})$ on that bracket is the stored ball.
Before any entry was written, the values for $3\leq k\leq7$ were compared with the ten-digit table of Dietzfelbinger, Goerdt, Mitzenmacher, Montanari, Pagh and Rink, and every value for $3\leq k\leq12$ was compared with the independent $2$-core characterisation in (3).