@Mac Dude 

Yes, it would be nice to be ensured that a+b-b+c  is not going to be interpreted as a+c+c.

@Christopher2222 

I was told that the bug exists in Maple 11 (classic).

@Christopher2222 

Unfortunalely the bug exists in all versions.

The bug seems to be very old; I tested in Maple 14 (2010) and it is present!

@Carl Love 

Since f^(1/p) is needed only when f' = 0 ==> all the exponents are multiples of p, so their division by p is ok.

And  r = c^(1/p) is simple when p=3: r=0, 1 or 2; actually r=c (by Fermat).

It is not clear: how many equations do you have?
As I understand, there are n = i+1 unknowns: S[1],...,S[i+1].
Do you use the letter i  only for n-1? You should restate mathematically the problem.

@Axel Vogt 

The beauty of such solutions largely compensates the absence of a direct answer from Maple.
And that's why mathematicians must exist!

@Markiyan Hirnyk 

It is interesting to note that with OP's solve, Maple does itself the conversion to rationals but after that it enters some huge computations with lots of memory and Windows freezes.

@Josolumoh 

I just said that all the moments (including the mean M[1]) are infinite:

M[k] = infinity, so, a "very" closed form!

 Manipulating the integral it is possible to compute this in Maple.
 Unfortunately I encountered some other bugs but they are probably related to the first one. 

Note also that Maple finds a correct (continuous) antiderivative but is not able to compute the limit at infinity.

 Text version:

(1/4)*sqrt(2)*sqrt(1-I)*(1+I)^(3/2)*EllipticE((1/2)*sqrt(2)*sqrt(1+I))+(1/4+1/4*I)*sqrt(2)*sqrt(1-I)*EllipticE((1/2)*sqrt(2)*sqrt(1+I), (1/2)*sqrt(2)-(1/2*I)*sqrt(2))+(-1/2-1/2*I)*EllipticE((1/2)*sqrt(2)-(1/2*I)*sqrt(2), (1/2)*sqrt(2)*sqrt(1+I))-(1/4*I)*sqrt(2)*sqrt(1-I)*EllipticF((1/2)*sqrt(2)*sqrt(1+I), (1/2)*sqrt(2)-(1/2*I)*sqrt(2))+(1/4+1/4*I)*EllipticF((1/2)*sqrt(2)-(1/2*I)*sqrt(2), (1/2)*sqrt(2)*sqrt(1+I))+1/2-(1/8)*sqrt(2)*sqrt(1-I)*(1+I)^(3/2)*EllipticK((1/2)*sqrt(2)*sqrt(1+I));

@tomleslie

It's worth mentioning  that for b[i]=0, the equations are satisfied by any P.

@Markiyan Hirnyk 

All the displayed digits in my comment are correct. It is easy to verify by solving dy/dx=0 with high precision.

@Markiyan Hirnyk 

Actually, with your displayed precision it is

2.0005898444608  at  x = 1.0003097363981.

OK, I have missunderstood the problem. I did not look for the matrix A. 

Actually A is not essential for the content; the author could have started with A without mentioning the method to produce it.

@Markiyan Hirnyk 

Yes. Proof: