Consider the following code:
c0 := -1.24; b := 0.125; hh[{x_, y_}] := {x^2 + c0 - b*y, x}; hh2[{x_, y_}] := hh[hh[{x, y}]]; dhh := Transpose[{D[hh[{x, y}], x], D[hh[{x, y}], y]}]; dhh2 := Transpose[{D[hh[hh[{x, y}]], x], D[hh[hh[{x, y}]], y]}]; fixpoints = NSolve[hh[{x, y}] == {x, y}, {x, y}]; period2cycle = NSolve[hh[hh[{x, y}]] == {x, y}, {x, y}]; center = {x, y} /. fixpoints[[1]]; centerOther = {x, y} /. fixpoints[[2]]; pushvector = Eigensystem[dhh /. fixpoints[[1]]][[2,1]]; pushvectorOther = Eigensystem[dhh /. fixpoints[[2]]][[2,1]]; sink = {x, y} /. period2cycle[[1]]; pushvector2 = Eigensystem[dhh /. fixpoints[[2]]][[2,1]]; maxIter = 63; maxXXSize = 600; maxYYSize = 600; smalldelta = 0.05/maxXXSize; newfun[{{x_, y_}, n_}] := {hh[{x, y}], n + 1}; newtest[{{x_, y_}, n_}] := Abs[x] < 10000 && Norm[{x, y} - sink] > 0.001; newcolorassign[{{x_, y_}, n_}] := Which[Abs[x] < 3, {0, 0, 0}, Mod[Floor[Arg[x]/2^16], 2] > 0, {0.5, 0.5, 0.5}, True, {1, 1, 1}]; With this setup one can execute the following:
Graphics[{Raster[Table[newcolorassign[NestWhile[newfun, {center - pushvector*smalldelta*I*(k - maxXXSize/2 + I*(j - maxYYSize/2)), 0}, newtest, 1, maxIter]], {k, 1, maxXXSize}, {j, 1, maxYYSize}]], PointSize[12/maxXXSize], White, Point[{(maxXXSize - 1)/2, (maxYYSize - 1)/2}], PointSize[4/maxXXSize], Black, Point[{(maxXXSize - 1)/2, (maxYYSize - 1)/2}]}] In a notebook this will produce a picture but because you specify the parameters at the beginning, to avoid confusion you should start a new notebook for each parameter pair (b,c0). Now to get around this problem I tried to generalise the code to the following.
hh1[{x_, y_}, b_, c_] := {x^2 + c - b*y, x}; hhhh[{x_, y_}, b_, c_] := hh1[hh1[{x, y}, b, c], b, c]; dhhh[b_, c_] := Transpose[{D[hh1[{x, y}, b, c], x], D[hh1[{x, y}, b, c], y]}]; dhhh2[b_, c_] := Transpose[{D[hhhh[{x, y}, b, c], x], D[hhhh[{x, y}, b, c], y]}]; fixpts[b_, c_] := NSolve[hh1[{x, y}, b, c] == {x, y}, {x, y}] per2cycle[b_, c_] := NSolve[hhhh[{x, y}, b, c] == {x, y}, {x, y}] centre[b_, c_] := {x, y} /. fixpts[b, c][[1]]; centreOther[b_, c_] := {x, y} /. fixpts[b, c][[2]]; pushvect[b_, c_] := Eigensystem[dhhh[b, c] /. fixpts[b, c][[1]]][[2,1]]; pushvectOther[b_, c_] := Eigensystem[dhhh[b, c] /. fixpts[b, c][[2]]][[2,1]]; snk[b_, c_] := {x, y} /. per2cycle[b, c][[1]]; pushvect2[b_, c_] := Eigensystem[dhhh[b, c] /. fixpts[b, c][[2]]][[2,1]]; maxIter = 63; maxXXSize = 600; maxYYSize = 600; smalldelta = 0.05/maxXXSize; newfun[{{x_, y_}, n_}, b_, c_] := {hh1[{x, y}, b, c], n + 1}; newtest[{{x_, y_}, n_}, b_, c_] := Abs[x] < 10000 && Norm[{x, y} - snk[b, c]] > 0.001; newcolorassign[{{x_, y_}, n_}] := Which[Abs[x] < 3, {0, 0, 0}, Mod[Floor[Arg[x]/2^16], 2] > 0, {0.5, 0.5, 0.5}, True, {1, 1, 1}]; unstablepic[b_, c_] := Graphics[ {Raster[Table[newcolorassign[NestWhile[newfun[#1, b, c] & , {centre[b, c] - pushvect[b, c]*smalldelta*I* (k - maxXXSize/2 + I*(j - maxYYSize/2)), 0}, newtest[#1, b, c] & , 1, maxIter]], {k, 1, maxXXSize}, {j, 1, maxYYSize}]], PointSize[12/maxXXSize], White, Point[{(maxXXSize - 1)/2, (maxYYSize - 1)/2}], PointSize[4/maxXXSize], Black, Point[{(maxXXSize - 1)/2, (maxYYSize - 1)/2}]}]; Now unstablepic[0.125,-1.24] produces the same output as the original code but is significantly slower. To see this more quickly, one can take maxXXSize = 60 and maxYYSize = 60 instead of 600. I'm not sure why the generalised code is taking so much longer than the original. Is there a way to make this quicker?
