@David Sycamore I found a counterexample to my Conjecture 2 as stated. I suspect that these counterexamples are extremely rare, and perhaps there's only a finite number of them; and I'll restate the Conjecture in a form that I think is very likely true (as the Riemann hypothesis is very likely true).

I think that the smallest counterexample is 81345, whose prime factor set is P:= {3, 5, 11, 17, 29}. It's easy to prove that this has exactly one solution (so not an infinite number): P has a complete set of remaiders mod 5. In order, they are {3, 0, 1, 2, 4}. So, for any e > 0, P +~ e necessarily contains a multiple of 5, which may be 5 itself. Indeed, using e=2 yields {5, 7, 13, 19, 31}, which are all prime, so that's a solution. But using any other e, you'll get one member that's a multiple of 5 other than 5 itself, so it won't be prime.

Note that I carefully chose the primes in P to not have a complete set of remainders mod 3.

Here's the restated conjecture: Conjecture 2: Let P be a finite set of primes which does not have a complete set of remainders with respect to any of its members. Then there are an infinite number of positive integers e such that P +~ e contains only primes.

This modified conjecture is still a generalization of the famous twin-prime conjecture.

Here are some conjectures related to this:

Conjecture 1: Every case where there is no solution can be verified by the method that I showed.

Conjecture 2: In every case where there is a solution, there are an infinite number of solutions.

The second conjecture is a generalization of the famous twin-prime conjecture; so, it may be a consequence of the Riemann hypothesis. 

I applied my proof technique to all 465 no-solution cases in Tom's failure list. In all 465 cases, it is proven that there are indeed no solutions. This is intended to be run after generating Tom's list failure:

ProveNoSolution:= proc(n)
local P:= numtheory:-factorset(n), p:= min(P); 
   evalb(nops(irem~(P, p)) = p)
end proc
:  
(proven, unproven):= selectremove(ProveNoSolution, op~(1, failure));

MaplePrimes can't handle .maple workbook files. Can you convert your file to a .zip file and repost it?

What you've posted is not a differential equation.

@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.