@minhthien2016 Denote by S the max of the sum of the radii for n=14.

The sum of the radii in Maple is SMaple = 1.886... > SMma = 1.880..., so, SMaple is better.
It is easy to check the constraints and see that Maple configuration is valid, so: S >= SMaple > SMma.

Note that actually S >= 1.897 and Maple can easily obtain this.
I will probably post this later.

If you have an Answer but want it be labeled as  Comment, the readers could be confused. I'd prefer to see an answer be presented as such.

@sursumCorda  Is convert(..., simplifier=F)  documented?

@Axel Vogt  Only under some assumptions, e.g. real.

@aroche  If was chosen the branch  diff(y(x), x) = (x-2)^(2/3)  instead of diff(y(x), x) = ((x-2)^2)^(1/3).
But then y(x) is complex for x<2, which is probably not wanted.

@nm It is 2 only because you want to be so. For the exponent 2, Maple says 4 (as you have complained).

But note that y(x) + x = 0 ==> y'(x)=-1, so you have:

(x + y(x)) * (1+diff(y(x),x)) = 0  <==> (1+diff(y(x),x)) = 0 <==>  y(x) = -x+c (for some constant c). 

@nm Mathematica is not the supreme authority.

P.S. With your statement,  (x + y(x))^7 * (1+diff(y(x),x)) = 0  has 8 solutions.

@nm It's not in the standard form (or reducible to)  d/dx y(x) = F(x, y(x)).

The reason is obvious. Anyway, I see ODESteps as designed for standard simple ODEs. 

@nm No bug without Physics updates.

@minhthien2016 Insert f:=simplify(f) assuming real; before calling extrema. 

@Carl Love Very nice, vote up! Unfortunately the rose is not quite 3d: rotating the plot, several gaps are visible and the stem is just "separately added".

@Hullzie16  Adding the option numeric, dsolve produces a nrmerical solution,
Your system has a lot of floats, so, approximations. In such situations, an exact (analytic) solution is usually useless (even when it exists). 

@Ali Hassani 

seriestorec(s, A(n)) 
   ==> for a given series s, finds a "suspected" recurrence relation of its coefficients A(n)
rectodiffeq(rec, A(n), f(z))
   ==> for a recurrence rec, finds the ode satisfied by the "generating function"
           f(z) = Sum(A(n)*z^n, n = 0 .. infinity).
 

@dharr  When checking the minimal polynomial we must use a much higher precision.

Digits:=200:
evalf(expr) - fsolve(poly)[2]: evalf[15](%);
#                        1.37954446875327 e-38

(it should be near 1e-200).
So, the first poly is incorrect.

The second one (of degree 20) is ok.