I am exploring the quadratic map $x_{n+1} = x_n^2 + c$. If we plot the value of $x$ after many iterations on the vertical axis against the input $c$ value on the horizontal axis we get a classic bifurcation diagram. Going from right to left, the diagram starts with a single stable curve at $\frac{1}{4}$. At $-\frac{3}{4}$ it splits into two branches, and then at $-\frac{5}{4}$ it splits into four branches. This period-doubling happens again and again, over increasingly narrow intervals, until the plot dissolves into chaos a little beyond $-1.401$. The ratio of the width of one stable region to the next tends towards the Feigenbaum constant of approximately 4.669.
I am trying to find the exact values of these bifurcation points from period-4 onward, or at least a method of approximating them. Wikipedia gives approximate values of $c$ for some of the bifurcation points, to eight significant figures for periods 8 through 256 and thirteen significant figures for periods 512 and 1024, plus eleven significant figures for the limit as the period goes to infinity. It does not explain how these values are calculated, however, and I haven't been able to find that information elsewhere.
How do I calculate these bifurcation points myself? My math education extends through Calc 1; if this requires math beyond that, I'd appreciate as much detail as you can provide so I can follow along.
I know that the bifurcations happen for values of $c$ such that there exist points where $f^{2n}(x) - x$ and the first derivative of that expression are both zero. For example, $c = -0.75$ at the first bifurcation point, where $f^2(x) - x = x^4 + 2cx^2 - x + c + c^2$ and its derivative $x^3 + 4cx - 1$ both equal zero for $x = -0.5$. For the second bifurcation point at $c = -1.25$, $f^4(x) - x$ and the first derivative both equal zero for $x = \frac{-1 ± \sqrt{2}}{2}$.
I can't go far with that, though, because the function $f^{2n}(x) - x$ gets exponentially longer as $n$ increases, to the point that the software I have available refuses to work with $n ≥ 8$.
The contraction factor $f_c'(x_1)f_c'(x_1)···f_c'(x_{2^n})$ along the cycle $x_1,x_2,...,x_{2^n}$ changes, from right to left, from $1$ to $-1$. There at the point of the bifurcation the cycle doubles to $x_1,x_2,...,x_{2^n},x_1,x_2,...,x_{2^n}$, and the new cycle starts again with $(-1)^2=1$. Along the old path the $2^n$ cycle continues, but is now unstable. Close-by a stable $2^{n+1}$ cycle branches off, locally looking like a pitch-fork bifurcation.
So you get two equations to solve to determine the two variables $x,c$ to get to a bifurcation point, $$F_{2^n}(x,c)=f_c^{2^n}(x)=x$$ and $$\partial_xF_{2^n}(x,c)=\prod_{k=0}^{2^n-1} f_c'(f_c^k(x))=-1$$
This behavior can be explored by plotting the function $x\mapsto f_c^{n}(x)-x$ for several values of $c$ close to one line of fixed points. From the graph the 4-cycle around $c=-1.3$ has a fixed point close to $0$.
One can see that as the slope at the root in the 4-cycle map becomes steeper in falling direction, the graph in the 8-cycle map develops a cubic looking shape, which is typical for this type of bifurcation.
It now becomes helpful to compute the Taylor polynomial of at least 3rd degree $f_c^N(x_0+h)=a_0+a_1h+a_2h^2+a_3h^3+O(h^4)$ to represent this pattern in an accessible form.
For the above plots the coefficients at the fixed points of the $4$-cycle are $$\small\begin{array}{c|l|l} c&\text{4-cycle coefficients}&\text{8-cycle coefficients}\\ \hline -1.300&[ 0.01943, 0.18054, 4.62595, -1.02845]&[ 0.01943, 0.03260, 0.98597, 7.53528]\\ -1.325&[-0.02352, -0.24406, 5.15920, 1.20675]&[-0.02352, 0.05957, -0.95184, -13.30452]\\ -1.350&[-0.05994, -0.67882, 5.48437, 3.01329]&[-0.05994, 0.46079, -1.19572, -43.82344]\\ -1.375&[-0.09219, -1.12390, 5.69108, 4.59867]&[-0.09219, 1.26316, 0.79252, -84.49985]\\ -1.400&[-0.12145, -1.57948, 5.82748, 6.09527]&[-0.12145, 2.49476, 5.33379, -140.92247]\\ \end{array}$$
Simply applying the root-finding idea works for the initial doubling cycles up to the $32\to64$ transition. After that the remaining space becomes so narrow and the iteration map too oscillating that strange defects in my implementation lead to paradox results, for instance graphing the leading coefficient function shows graphs that do not take negative values.
$$\small\begin{array}{c|l|l} N & \text{bifurcation parameter $c$ for $N\to2N$ cycles}\\ \hline 1 & -0.750000000000 \\ 2 & -1.250000000000 \\ 4 & -1.368098939391 \\ 8 & -1.394046156601 \\ 16 & -1.399631238874 \\ 32 & -1.400828742371 \\ \end{array}$$