@Ronan You have a more subtle mistake, corrected but not mentioned by the responders, such that just replacing 0 with NULL is not enough. Try to find and explain it!

This reminds me that for a polynomial with random coefficients, its roots tend to lie near the unit circle.

p:=randpoly(x, degree=200, dense):
plots:-complexplot([fsolve(p, complex)], style=point, scaling=constrained)

Maybe LaTeX should accept empty matrices (even if they do not make much sense).

@Tokoro For n=5, R is the same R=8/sqrt(15) and the edges of the 5-gon: 2,3,2,3,2.

@Tokoro  I had in mind the existence of any number of points at integer distances on a circle or sphere (not a trivial result). The minimal radius depends of course on the dimension. But you could have formulated the problem in R^d in the first post.

@nm You can simply define then

Dy := t -> D(y)(t);
D2y := t -> (D@@2)(y)(t);

Another approach would be to write a proc, converting (via readline) ode.mpl to trueode.mpl, by replacing y'(u) with D(y)(u) etc.

Why don't you add this version to the previous very similar one? Probably a more strict admin will delete this.
You have received answers to the previous question, but you seem to ignore them.
Note that the answer to the circle variant is automatically an answer to the sphere one!

@Carl Love The distance between any pair of points should be integer!

@bstuan t = 3 is not an assignment; it is an equation and the variable t remains free.
The assignment operator is :=.
So, only after t := 7;  the exprassion 10*t+5; will be evaluated to 75.

@tomleslie The situation is even worse because both solutions are wrong; they are valid only for b<=0.

restart;
sys:=[a + d = 0, -c + d = 0, a - b = sqrt(a^2 + b^2 + c^2)]:
sol:=solve(%) assuming real:
eval(sys,sol[1]), eval(sys,sol[2]) ;

        [0 = 0, 0 = 0, -b = (b^2)^(1/2)], [0 = 0, 0 = 0, -3*b = 9^(1/2)*(b^2)^(1/2)]

Download parab-rat.mw

@Preben Alsholm I mean: cannot simplify to  2*Pi/(p*sin(Pi/p))

Edit. Converting first to GAMMA it works.

@Preben Alsholm A human would reduce the integral to a Beta function using the change of variables x^p = t.
Unfortunately, Maple cannot compute it. 
Actually, including the assumption p>1 (mandatory for convergence), the verbose int finds this:

[cook = 2*Beta(1/p, -1/p + 1)/p, FAILS = (distribution, piecewise, series, o, polynomial, ln, lookup, ratpoly, elliptic, elliptictrig, meijergspecial, improper, asymptotic, ftoc, contour), ftocms = 2*GAMMA((p + 1)/p)*GAMMA((p + 1)/p - 2/p), meijerg = 2*Pi/(p*sin(Pi/p))]
 

@Carl Love Thank you. I was aware that the proc can be improved but I have chosen the simplicity, and anyway a Dynkin system is used mainly for infinite families and in this case Maple is useless.
My opinion is that a Maple code implementing a more complicated algorithm (not the case here) should be presented in two flavours: a simple one and an optimized (speed and/or memory) version.

@JAMET But why?