I'm beginner and do not how to plot the Basins of attraction by the methods

Question

For a complex polynomial $P(z)$, define $T$ by $T(z)=z-\frac{P(z)}{P^{\prime}(z)}$, where $P^{\prime}(z)$ is the first derivative of $P(z)$. We consider Newton method for finding the roots of polynomial $P(z)$, which is given by the recursive formula

z_{n+1}=T(z_{n}), \ n\in\mathbb{N},\ \ \ \ (16)

where $T$ is defined as above. Similarly, we set a new method (and called it S iteration) where $T$ is defined as above, by

y_{n}=(1-\beta_{n})z_{n}+\beta_{n}Tz_{n}
z_{n+1}=(1-\alpha_{n})Tz_{n}+\alpha_{n}Ty_{n}, n\in\mathbb{N}, (17)

Now If the sequence $\{z_n\}_n\in\mathbb{N}$ (the orbit of the point $z_1$) converges to a root $z^{\ast}$ of the polynomial $P$, then we say that $z_1$ is attracted by $z^{\ast}$. The attraction basin of the root $z^{\ast}$ of the polynomial $P$ is the set of all starting points $z_1$ which are attracted by $z^{\ast}$. If the test polynomial is

P(z) = z^5 + z^4 + z^3 + z^2 + z,  (18)

For this test polynomials $P(z)$ above, we consider square domain is centered at the origin in the complex plane as follows

D=[-10,10]\times[-10,10].

If $\alpha_n=\beta_{n}= 0.8$ for all $n \in\mathbb{N}$ in S iteration method (17), To generate the basins of attraction and to study the dynamics of methods (16) and (17), we divide the domains $D$ into 250 × 250 grids. If the sequence $\{z_n\}_n\in\mathbb{N}$ attempts a root of polynomial $P$ with accuracy of $10^−4$ in number of iterations $n \leq 13$, then the converging point $z_0$ is colored in a color assigned to $n$; otherwise, the point is colored in white. Figure 1a, b (see the image in attach or in the link for clear picture https://link.springer.com/article/10.1007/s00500-020-05038-9) represents the obtained basins of attraction by Newton method (16) and S-iteration method (17) for the polynomial $P$ defined in (18).

Now the attached pictures are generated by matlab and I want to generate the same pictures in mathematica.

Answer 1

For Newton's method, this could look as follows:

p = z^5 + z^4 + z^3 + z^2 + z
znew = Together[z - p/D[p, z]]

cf = With[{
    p = p,
    num = N@HornerForm[Numerator[znew]],
    denom = N@HornerForm[Denominator[znew]]
    },
   Compile[{{z0, _Complex}, {maxiter, _Integer}, {\[Epsilon], _Real}},
    Block[{iter, z},
     z = z0;
     iter = 0;
     While[Abs[z] > \[Epsilon] && iter < maxiter,
      iter++;
      z = num/denom
      ];
     iter
     ],
    CompilationTarget -> "C",
    RuntimeAttributes -> {Listable},
    Parallelization -> True
    ]
   ];

n = 1000;
gridpts = Tuples[Subdivide[-10., 10., n - 1], 2] . {1., 1. I};
result = cf[gridpts, 13, 1. 10^-4];

ArrayPlot[Partition[N[result], n],
 ColorFunction -> "SunsetColors"
 ]