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.

So, the problem reduces to show that
F <= 2   for  0<z<1
where
F:=(1-z)^(2*sqrt(z+1))+(z+1)^(2*sqrt(1-z))

Note first that the Maple solution (based on numerical methods)  is not acceptable from a mathematical point of view.

F is analytic in (-1,1), and also in the unit disk of C.
It has the power series

  

It would be enough to prove that all the coefficients c_k  of z^(2k)  are >0, (k>0).
So, the problem is reduced from "continuous" to "discrete" but still not easy.
It is obvious that F<2 in (0,a)  for a>0 small enough.
It is also easy to see that F<2 for a >= a0 =  e^(LambertW(4 (e)^2)-2)-1 = 0.613... = the root of  the derivative of the second term of F.

So, it remains the interval (0,a0).  This could (?) be treated by trying to obtain some estimates for the coefficients c_k of the power series.

Edit. The mentioned method works. Here are the details:

a) Using the Cauchy estimates,   |c_k| < 3, for any k.

b) Using a few terms of the series, for 0 < z <= a = a0:

f(z) < 2-(2/3)*z^4-(9/160)*z^6-(1/126)*z^8 + 3 * z^10*(1 + a^2 + a^4 + ...)
<= 2-(2/3)*z^4-(9/160)*z^6-(1/126)*z^8 + 3*z^10/(1-a^2)

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

c) The polynomial with rational coefficients p  is <= 2  for z in [0,a] subset [0,65/100].

@John Fredsted 

OP says 2's complement but means 1's ?

You should use/ask Maple code only after there is a clear mathematical formulation.

So, you have a curve given implicitely f(u,v)=0 [after your change of coordinates] and a point (u0,v0) on it, f(u0,v0)=0.

Do you want an approximate parametric representation of the curve near (u0,v0)?

u = u0 + a1 z + a2 z^2 + ...
v = v0 + b1 z + b2 z^2 + ...

Why don't you start with a parametric representation? Your initial curve is rational, so Maple can find such a representation if you don't know one.
Is this your problem?

@Carl Love 

I was not very interested in speed or length; it was simply a joke, knowing the "compactness" of your code (which I like, BTW).

But I do know that f(g(x)) is usually faster than (f@g)(x)

( see showstat(`@`))

@Carl Love 

I am puzzled by the fact that you did not use

(parse@cat)(L[]);    # one character shorter

or
parse(cat(L[]));      # two!

(also faster, because it avoids @).

It is not in your style!   :-)

@Markiyan Hirnyk 

It's not wrong at all. Verify!

What could be wrong in my 1 line posted solution:

g(c,d) >= max( f(1), f(-1)) = max (-1-d, 1+d) = abs(1+d) >= 0

?

@Markiyan Hirnyk 

If you can't understand the convexity stuff, just omit it, it is only for the [0,1] case.