Feeling bored and restless. Dug up some old code I had started to write a while ago when I thought that eigenvector methods would be a great way to render the feigenbaum diagram. Actually, it doesn't really work that well at all:



Any bigger and it starts yelling "ComplexInfinity" at me. Stupid mathematica.

ig = Function[{a,p}, Re[1/2 - Sqrt[1/4 - p / a]]];

d = Function[{a,n},
arr = Table[0, {i,1,n}, {j,1,n}];
For[i=0.0, i < n, i++, 
 y0 = i / n; 
 y1 = (i+1) / n;
 x0 = ig[a,y0];
 If[y1 / a > 0.25,
   Block[{}, (* If we're about to leave the range of the function *)
    xi = Floor[x0 * n];
    v = (xi + 1) - x0 * n;
    If [xi < n / 2, arr[[xi + 1,i + 1]] = v];
    For[xj = xi + 1, xj < n/2, xj++, arr[[xj + 1, i + 1]] = 1]; 
    Break[];
   ],
   Block[{}, (* If not *)
    x1 = ig[a,y1];
(*    Print[x1]; *)
    xi = Floor[x0 * n];
    While[x0 < x1,
(*      Print["---> ", x0," ", x1, " ", xi];*)
     If[(xi + 1) / n > x1, v = (x1 - x0) * n, v = (xi + 1) - x0 * n];
     arr[[xi + 1,i + 1]] = v;
     xi++;
     x0 = xi / n;
    ]
   ]
 ]
];
For[i = 0, i <= n / 2, i++, arr[[n - i]] = arr[[i + 1]]];
Eigenvectors[Transpose[arr],1]
];

make = Function[{xn,yn},
Block[{z,w},
data = {};
For[z = 0, z < xn, z++,
val = 3 * z / xn + 1;
ww = (d[val, yn])[[1]];
ww = ww / Total[ww];
data = Append[data, ww];
];
img =  (Floor[255 - 255 * #])& /@ Reverse[Transpose[data]];
(* Export["b.dat", Reverse[Transpose[data]], "Table"]; *)
Export["a.pgm", Join[{"P2", xn, yn, 255}, Flatten[img]],"List"];
]];