Conclusion
The special point is at $$c = 0.2294487977025897\ldots$$ satisfying $$f_c(f_c(z_0(c))) = f_c(z_0(c))$$ where $$z_0(c) = \sqrt[3]{-\frac{2}{5c}}$$ is the non-zero real critical point. Its filled-in Julia set has a single component.
Finding the point
Script for (wx)Maxima computer algebra system:
f(c, z) := z^2 + c * z^5;
critical_points : solve(diff(f(c, z), z) = 0, z);
front(l) := firstn(l, length(l) - 1);
solutions : unique(xreduce(append, front(map
( lambda([cp], solve([at(f(c, f(c, z)) = f(c, z), cp)], [c]))
, critical_points
))));
find_root(solutions[1], c, 0.2, 0.3);
Being more explicit:
solve([at(f(c, f(c, z)) = f(c, z), critical_points[3])], [c]);
%[2];
% * 9375 * 2^(1/3);
% + 81*2^(4/3);
%^2;
expand(%);
at(%, [c = x^3]);
solve(%, x);
%[2];
at(%, [x = c^(1/3)]);
gives the equation
$$5^{35/3}c^{5/3}-5^{10}\cdot3\cdot2^{2/3}c-5^5\cdot3^4\cdot2^{8/3} = 0$$
whose real root $c \in [0.2, 0.3]$ seems to match the solution found in the previous code block.
So, yes, you need to solve a quintic (a polynomial in $x = c^{1/3}$).
Finding the number of components
Consider the long-term behaviour of (one of) the other non-zero complex conjugate pair of critical points to know how many components the Julia set has for this specific $c$. Using SageMath I found that the other non-zero critical point $z_1(c)$ is pre-periodic, with
$$f_c(z_1(c)) = f_c(f_c(f_c(z_1(c))))$$
so it remains bounded. As no critical points escape to infinity, the Julia set for this $c$ value has a single component.
Script for SageMath:
x = polygen(QQbar)
f(c, z) = z^2 + c * z^5
five3 = QQbar.polynomial_root(x^3 - 5, CIF(RIF(0, 2), RIF(-1, 1)))
two3 = QQbar.polynomial_root(x^3 - 2, CIF(RIF(0, 2), RIF(-1, 1)))
P = five3^35 * x^5 - (5^10 * 3 * two3^2) * x^3 - 5^5 * 3^4 * two3^8
c0 = QQbar.polynomial_root(P, CIF(RIF(0.6, 0.7), RIF(-0.1, 0.1)))^3
print(c0)
z0 = QQbar(- (2 / (5 * c0))^(1/3))
z1 = QQbar((- 2 / (5 * c0))^(1/3))
print(z0, f(c0, z0), f(c0, f(c0, z0)))
print(z1, f(c0, z1), f(c0, f(c0, z1)), f(c0, f(c0, f(c0, z1))))
print(all([f(c0, z0) == f(c0, f(c0, z0)), f(c0, z1) == f(c0, f(c0, f(c0, z1)))]))
The point is in the right place
From the left
Arbitrarily close nearby values of $c$ less than the special point give filled-in Julia sets with infinitely many components.
Considering small real perturbations in $c$ near the special point and non-zero real critical point, the difference in $z$ has the opposite sign to the difference in $c$. Now considering small real perturbations in $z$ near the eventually fixed point in the limit of the perturbation in $c$ going to $0$, if the perturbation in $z$ is initially positive, then it grows without bound, eventually escaping towards infinity.
Continuing with SageMath:
zs = f(c0, z0)
F(C, Z, c, z) = (f(C + c, Z + z) - f(C, Z)).expand()
print(F(c0, z0, c, 0)) # F has opposite sign to c
print(F(c0, zs, 0, z)) # F > 2z if z > 0
From the right
It remains to be shown that arbitrarily close nearby real values of $c$ greater than the special point give filled-in Julia sets with a single component.
If we have $|z_{n+1}| < |z_n| < \epsilon$, then iterations are certainly bounded as $n \to \infty$. Substituting, we get
$$|\epsilon^2 + c \epsilon^5| < \epsilon^2 + |c| \epsilon^5 < \epsilon $$
And if $|c| > \frac{2}{5\epsilon^3}$ then $z_0 < \epsilon$. So
$$\frac{2}{5\epsilon^3} < |c| < \frac{1 - \epsilon^3}{\epsilon^4}$$
whence $$5 \epsilon^3 + 2 \epsilon < 5,$$ that is $$\epsilon < \epsilon_{\text{max}} \approx 0.86755924139011\ldots$$
and
$$|c| > R_{\text{min}} \approx 0.6125796570109\ldots$$
That still leaves a big gap of positive real $c$ that haven't been proven to remain bounded.
Visualisation
Here's the "Mandelbrot set" plot of the $c$ plane with the proven-bounded part in red, points that definitely escaped to infinity in white (these have filled-in Julia sets with infinitely many components, the rest have a single component), and the remainder in black:
Here's a zoom-in of the previous image with a grid overlaid in green (spacing 0.1 units, the axes are slightly thicker); the special $c$ is at the right-most corner of the white region:
Here's an animation of the filled in Julia set (black) for the special $c \pm \frac{1}{10}$.
Parameter plane GLSL fragment shader source code snippet for FragM (fork of Fragmentarium):
#version 330 core
#include "TwoD.frag"
const float pi = 3.141592653;
vec2 cMul(vec2 a, vec2 b)
{
return vec2(a.x * b.x - a.y * b.y, a.x * b.y + a.y * b.x);
}
vec2 cConj(vec2 a)
{
return vec2(a.x, -a.y);
}
vec2 cDiv(vec2 a, vec2 b)
{
return cMul(a, cConj(b) / dot(b, b));
}
vec2 cExp(vec2 a)
{
return exp(a.x) * vec2(cos(a.y), sin(a.y));
}
vec2 cLog(vec2 a)
{
return vec2(log(length(a)), atan(a.y, a.x));
}
vec2 cCbrt(vec2 a, int n)
{
return cExp((cLog(a) + vec2(0, float(n) * 2.0 * pi)) / 3.0);
}
vec3 color(vec2 c, vec2 dx, vec2 dy) {
bvec3 escaped = bvec3(false);
for (int j = 0; j < 3; ++j)
{
vec2 z = cCbrt(cDiv(vec2(-2.0/5.0, 0.0), c), j);
for (int i = 0; i < 1000; ++i)
{
escaped[j] = escaped[j] || !(dot(z, z) < 1.0e6);
if (escaped[j]) break;
vec2 z2 = cMul(z, z);
vec2 z4 = cMul(z2, z2);
vec2 cz4 = cMul(c, z4);
z = cMul(z, z + cz4);
}
}
vec2 g = c;
if (abs(g.y) < dy.y) return vec3(0.0, 1.0, 0.0);
if (abs(g.x) < dy.y) return vec3(0.0, 1.0, 0.0);
g *= 10.0;
g += 0.5;
g -= floor(g);
g -= 0.5;
if (abs(g.y) < 2 * dy.y) return vec3(0.0, 1.0, 0.0);
if (abs(g.x) < 2 * dy.y) return vec3(0.0, 1.0, 0.0);
return mix(vec3(escaped), vec3(1.0, 0.0, 0.0), float((length(c) > 0.61257965701099)));
}
Filled-in Julia set GLSL fragment shader source code snippet for FragM:
#version 330 core
#include "TwoD.frag"
const float pi = 3.141592653;
vec2 cMul(vec2 a, vec2 b)
{
return vec2(a.x * b.x - a.y * b.y, a.x * b.y + a.y * b.x);
}
vec3 color(vec2 z, vec2 dx, vec2 dy) {
vec2 c = vec2(0.2294487977025897 - 0.1 * cos(2.0 * pi * time), 0.0);
bool escaped = false;
for (int i = 0; i < 1000; ++i)
{
escaped = escaped || !(dot(z, z) < 1.0e6);
if (escaped) break;
vec2 z2 = cMul(z, z);
vec2 z4 = cMul(z2, z2);
vec2 cz4 = cMul(c, z4);
z = cMul(z, z + cz4);
}
return vec3(escaped);
}
Reproducibility
Try the SageMath scripts online here. Today the output is:
0.2294487977025897? + 0.?e-37*I
-1.203533181541392? 0.8690952714426872? 0.8690952714426872?
0.6017665907706960? + 1.042290309512355?*I -0.4345476357213436? + 0.752658583378300?*I -0.4345476357213436? - 0.752658583378300?*I -0.4345476357213436? + 0.752658583378300?*I
True
-2.525168374948163?*c
0.2294487977025897?*z^5 + 0.9970643256076521?*z^4 + 1.733087781419605?*z^3 + 2.506218395826876?*z^2 + 2.392714185671939?*z
I used the Debian package for maxima, version 5.46.0-8 (not the latest version available in the Bookworm repositories).
I used https://sagecell.sagemath.org as interface to SageMath.
Images were rendered with Fragmentarium 2.5.7.221026 Digilantism with my example files (include path: Examples/Include;Examples/Claude).
Best Answer
You wrote:
Sure! This little web app generates images of Julia sets for functions of the form $f_c(z) = z^2+c$ using the modified inverse iteration algorithm as described in "The Beauty of Fractals" and "The Science of Fractal Images". Viewing the page source, it's pretty easy to find this Javascript code for the program.
The basic ideas behind inverse iteration, it's modification appear, and implementation for Mathematica were described in this 1998 paper. That program has been been built into the Mathematica kernel since V10. One thing that's nice about that version is that (with all the algebraic machinery that Mathematica has to offer), it's not hard to write a program that works for higher order polynomials or rational functions.
Another technique to get the boundary is called "boundary scanning". This is particularly useful for certain hard to generate Julia sets, as described in this answer.
Here's a comparison of those two algorithms for $z^2 - 0.77967939051932 + 0.11124251677182495i$:
Inverse iteration
Boundary scanning