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"]; ]];