@Markiyan Hirnyk 

I know, but I do not see anything wrong in my previous comment. On the contrary, the information should be useful for anyone using extrema or similar commands.

@Markiyan Hirnyk 

In this example the situation is the same. extrema returns the values (min and max) for the functions at the critical points (but complex included). In this case a single expression but it is multi-valued being a RootOf.
The user must select only the real critical points and then take the min/max of the corresponding values.

Note that the general situation is complicated even for analytic functions of a single variable because usually (including this example) the set of critical points is infinite; allvalues gives only a finite set.
RootFinding[Analytic]  finds easily many other critical points!

Edit. That's why maths is always important when using any CAS!

@Markiyan Hirnyk 

Yes, but this happens because for one of the complex critical points, the value of the function is negative, so < 9/4.

@Markiyan Hirnyk 

But {r = 1, s = 1}  is visible in your worksheet!

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