One idea for computing the Feigenbaum constant is as follows: $\delta$ is the unique expanding eigenvalue for period-doubling renormalization at its fixed point $F$. So take any real-analytic 1-parameter family of univalent maps that is transverse to the stable manifold of renormalization (for example, the real quadratic family). Then the iterated renormalizations of this 1-parameter family will converge to the unstable manifold of renormalization. The unstable manifold is of course mapped to itself by renormalization and the derivative at the fixed point will be $\delta$.
It should be possible to implement this numerically, by keeping track of the power series for a real-analytic family of real-analytic univalent map, and applying the renormalization by replacing $f_\lambda$ with $f_\lambda \circ f_\lambda$ (and suitably rescaling). Because renormalization acts as a contraction on these 1-parameter families, this procedure should be computationally stable, and provide a number of digits proportional to the number of times one renormalizes.
One classical reference for the stable/unstable manifold picture for period doubling renormalization is Iterated Maps on the Interval as Dynamical Systems by Collet and Eckmann.
mpfun
). $\endgroup$