@JAMET These type of questions are better  to be presented mathematically and followed by the worksheet with your approach. Otherwise, the math can be skewed by errors in your code and it takes more time for someone to help you. Many people here will not waste their time just because you refuse a proper presentation of the problem.

@sand15  I do not see what is the numerical benefit of using r[i] >= - eps, with eps =1e-8.
It is obvious that mathematically the conditions r[i] >=0  or  r[i] >= +eps  are equivalent (i.e.  S is exactly the same).

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