@Markiyan Hirnyk 

The issue here is of a different nature: extrema computed the critical points but wrongly included also the complex ones.

minimize/maximize seem to be more reliable commands for unconstrained min/max. It works for:

f:=(x*y+x*z+y*z)*(1/(x+y)^2+1/(y+z)^2+1/(z+x)^2);
g:=eval(f,[x=1,y = s^2, z = r^2]);
minimize(g,location);

@Klaroline 

Still not clear.
If you want the values of a satisfying g(c)<1 for each critical point c of g
then you get the empty set for the original function g(x)=exp(x^2*(x-a)).

You should start by explaining the exact meaning for range of "a".

@Klaroline 

@maniaerox 

That is exactly what a parametric plot is doing. Try it with your expressions.

You could add options such as:
plot( [E,T, w=0..2*Pi],labels=["E","T"], view=[-1/2 .. 2, -2..2]);

The RootOfs should come from the critical points.

f:=(a+b+c+d)^2*(a*b+a*c+a*d+b*c+b*d+c*d)^2-144*(a^2+b^2+c^2+d^2)*a*b*c*d:
g:=expand(eval(f,d=1)):
S:=[diff(g,a),diff(g,b),diff(g,c)]:
G:=Groebner[Basis](S,plex(a,b,c)):
p:=factor(G[1]);

fsolve(op(7,p));

fsolve(op(8,p));

@Preben Alsholm 
You have developped this interesting tool related to a recently posted system of ODEs.
Just curious: did you see in the past many such systems where unwanted expansions occur during dsolve?

@Axel Vogt 

He simply expressed the original symmetric function wrt other symmetric functions (those pij defined).  But he did not prove the inequality. I think that he suggested to look for other pij which could lead to success.

If we set at the beginning

Digits:=1000;

we obtain quickly the expected s:=0  (!!!).

@Markiyan Hirnyk 

1. There are minor inaccuracies in it: 1 = a < b < c < d should be replaced by a = 1, b >= a, c >= b, d >= c.

No, I said: we may suppose...  And we may, because once the inequality is proved for a<b<...  it holds in general (by continuity). 

2. It was a real pleasure for me to discuss with you without attacks...
Same for me, I hope that you will maintain this direction. 

PS. The code can be executed immediately by copy&paste, so no need for a .mw.

@Markiyan Hirnyk 

It works in Maple 2015. Your worksheet contains only the error message, not the input.

If you do not have an input error, I will rewrite the select.

Edit. I edited simplify to expand.
Probably Maple 2016 simplifies differently. Is it ok?

for the inequality

(a,b,c,d>0).

The proof is mathematically acceptable because it can be verified by hand.

The function being symmetric and homogeneous we may suppose 1 = a < b < c < d.
f:=(a+b+c+d)^2*(a*b+a*c+a*d+b*c+b*d+c*d)^2-144*(a^2+b^2+c^2+d^2)*a*b*c*d:
g:=expand(  eval(f, [a=1,b=1+x, c=1+x+y, d=1+x+y+z])  ):  # x,y,z>0
g1:=select( u -> coeffs(u)>0, g);

g2:=select( u -> coeffs(u)<0, g);

The idea is to show that these negative terms can be combined with some terms in g1 such that their sum be positive.

A:= 60*x*z^3 + 36*z^3;
B:= 12*y*z^3+24*y*z^2;  #g2 = -(A+B)

AA:= 36*x^3*z^3+4*x^2*z^4+450*x^4*z+246*x^3*z^2+14*x^2*z^3+12*x*z^4+444*x^3*z+154*x^2*z^2+9*z^4+180*x^2*z+60*x*z^2+24*x*z+36*z^2;
AA-A; factor(%): p1:=op(2,%):
collect(p1,z);

p1 cannot have positive roots in z by Descartes' rule of signs.  Hence A <= AA.
##############################################
BB:= 4*y^6+12*y^5*z+13*y^4*z^2+6*y^3*z^3+y^2*z^4+64*y^5+160*y^4*z+140*y^3*z^2+50*y^2*z^3+6*y*z^4+112*y^4+224*y^3*z+100*y^2*z^2+96*y^3+144*y^2*z+48*y^2+48*y*z:
{op(AA)} intersect {op(BB)}, {op(AA),op(BB)} subset {op(g1)};

   
BB-B: factor(%): p2:=op(2,%):

collect(p2,y);

factor(coeff(%,y,0));

Hence p2 has nonegative coefficients --> p2>0 for y,z>0  ==> B <= BB
So,  g >= (AA-A) + (BB-B) >= 0.

Q.E.D.

@Axel Vogt 

I took M0 = sup{ |F(u)| : u in C, |u|<1 }. Then R=1  (actually, R--> 1, R<1).

M0 < M = 3.  (==> |c_k| < M)

To obtain M=3 some evaluations are needed. But it is possible to use a bigger bound increasing the number of terms.

For example one may use M0 < M = 10 (which is easy to obtain) and take

p = 2-(2/3)*z^4-(9/160)*z^6-(1/126)*z^8-(72013/2580480)*z^10 + 10  * z^12/(1-(7/10)^2);

Using your change of coordonates, your points on the (u,v) curve are all complex. Do you really want this?

precurve3 := 1-Q*beta+alpha*beta^3*(Q^2*beta^4-Q*beta^3-2*Q*beta^2+2*beta^2-beta+1)-alpha^2*beta^9*(1-beta):
sol := solve(eval([precurve3 = 0, diff(precurve3, beta) = 0], Q = 2.0), [alpha, beta]):
curve := subs([alpha = exp(u), beta = exp(v)], eval(precurve3, Q = 2)):
sol3 := [seq([u = ln(rhs(sol[i][1])), v = ln(rhs(sol[i][2]))], i = 1 .. 8)]:
s1 := solve(curve, u):
u0,v0:=eval([u,v],sol3[2])[]:
v=v0+z;

u=series(  eval(s1[2],v=v0+z), z);

I have anticipated your position.
I let (and invite) other mathematicians here to comment.
But I provided a solution in the edit of the comment "Just some ideas".

@Markiyan Hirnyk AFAIK, DS never finds exact solutions. It works purely numerical.
But I do not want to continue the discussion in this direction.