@patrickfitzpatrick Please let me know whether the all-solutions code that I posted works adequately for you. In particular, Does it produce all solutions in a reasonable time for all the cases that you're interested in?

@acer Unfortunately, workbook .maple files don't work on MaplePrimes.

@Carl Love The algorithm for the all-solutions problem and its Maple implementation are in an Answer below, and I'm quite pleased with them.

@patrickfitzpatrick Both Kitonum's and my code give one solution--not all solutions--for a given n. I am working on an algorithm to efficiently find all solutions. It should take me a few hours.  

@Masooma The version released in 2018 is called "Maple 2018". There was also a version released several years earlier (2014?) called "Maple 18". Which of these two do you have? The command kernelopts(version) will tell you.

Why do you have size=[600,300] if you want another size?

@Masooma The only part of your code above that makes sense is the first line:

F := proc (x) options operator, arrow; x-2*f(x)/(D(f))(x) end procyou

It's more commonly written simply as

F:= x-> x - 2*f(x)/D(f)(x)

All the stuff with algsubs does nothing, and I don't know what you're trying to accomplish. The thing that you wrote with a second derivative could be meaningful, but it doesn't seem relevant to your goal of iterating the first function F.

By the way, what version of Maple are you using? I could retrofit my code to an earlier version if you need that.

@Masooma In order for the order of convergence to be 2, m must be the known multiplicity of the root (many thanks to Acer and VV for providing clarity on this). So, if we use m=2, and the root has multiplicity 2, then the order of convergence will be 2. If you apply my order-of-convergence approximator (code given above) to such a situation, you'll get an order very close to 2. The polynomial examples that you gave earlier do not have roots of multiplicity 2:

(x-1)^3                                   m=3
(x^2+2*x+1)^5 = (x+10)^10   m=10

What size n are you interested in expressing as a sum of 4 squares? The algorithm that I gave below works reasonably efficiently (say, <= 16 milliseconds) for n of many hundreds of bits. It may get bogged down if n has two or more large, distinct prime factors on the order of, say, a hundred bits. This slowness is just due to finding the prime factorization of n. If you need to work with such n, I'll implement the Rabin & Shallit algorithm (from the paper mentioned below), which doesn't use the prime factorization of n and has average time complexity O(log(n)^2*log(log(n))).

Yes, it's easy. What do you want to be the time parameter of the animation--perhaps the standard deviation? You'd likely get better pedagogy---if that's your purpose---from an Explore command than from a dedicated animation.

@vv Yes, that's the paper that I was looking for. Thank you for posting it.

When I first read this Post, I quickly suspected that the order of the permutation would be a fairly easy closed-form function of its length and that the number of disjoint cycles would be small regardless of  length; however, neither of those seem to be true. Given that, it would seem odd that Order(P(n)) <= n for all n. Has someone proven that? Below, I've verified it up to n=4096.

Here are some explorations:

Download ArabicCipher.mw

I suppose the problem could be done in Maple. Matlab's fminsearch seems equivalent to parts of Maple's third-party DirectSearch package. I need details of "Pareto front" to proceed.

@Masooma I am happy to help further if you is still need it. Newton's method and its variations are an area of expertise of mine. I still don't see any merit to your method x-> x - 2*f(x)/D(f)(x) whether for ordinary or repeated roots.

It wouldn't let me attach the worksheet above. Here it is:
 

The histograms can be generated in the worksheet. There's some trouble with uploading them directly.

Download UniformRandPrime.mw