Optimization:-NLPSolve works very well and fast without the need of any simplification.

restart;
Digits:=15:
e:=(x,y) -> x^2 + 2*y^2 - 1:
h:=(x,y) -> (x-sin(x))^2 + (y-sin(y))^2 - 1:
constr:=
e(x1,y1)=0, h(x2,y2)=0, h(x3,y3)=0, r>=0,
D[1](e)(x1,y1)*(y1-y0) - D[2](e)(x1,y1)*(x1-x0) = 0,
D[1](h)(x2,y2)*(y2-y0) - D[2](h)(x2,y2)*(x2-x0) = 0,
D[1](h)(x3,y3)*(y3-y0) - D[2](h)(x3,y3)*(x3-x0) = 0,
(x1-x0)^2 + (y1-y0)^2 = r^2,
(x2-x0)^2 + (y2-y0)^2 = r^2,
(x3-x0)^2 + (y3-y0)^2 = r^2,
x3 >= x2 + 1/10:

sol:=Optimization:-Maximize(r, {constr},   x0=1..1.4, y0=1..1.4, x1=0..1, y1=0..1, x2=1..2, y2=1..2, x3=1..2, y3=1..2);

sol := [0.769863723979874459, [r = 0.769863723979874, x0 = 1.13805882552303, x1 = 0.699729170630851, x2 = 1.27446720511737, x3 = 1.89574141571275, y0 = 1.13805882552303, y1 = 0.505162888466810, y2 = 1.89574141571275, y3 = 1.27446720511737]]

p1:=plots:-implicitplot(e, -2..2, -2..2): p2:=plots:-implicitplot(h, -2..2, -2..2):
p3:=plottools:-circle(eval([x0,y0],sol[2]), eval(r,sol[2]), color=red):
p4:=plots:-pointplot(eval([[x0,y0],[x1,y1],[x2,y2],[x3,y3]],sol[2]), symbolsize=8, color=blue):
p5:=plots:-textplot(eval({seq}([x||i,y||i,P__||i], i=0..3), sol[2]), align={below, left}, color=blue):
plots:-display(p1,p2,p3,p4,p5);

data_new := DataFrame(data_2, rows=[1,2]);

For a system of m polynomial equations of (total) degrees d_1, ..., d_m in n unknowns over C (or other algebraically closed field), the standard result is given by Bezout' theorem:

For m=n, the number of solutions is either infinite or <= d_1* d_2 * ...* d_m.

For details see: Bézout's theorem - Wikipedia

restart;
# f := floor;
f := x -> x^2:
X := [seq(x, x=-3..3,0.02)]: n:=nops(X):
Y := map(f, X):
C := Array(1..n, i -> X[i], datatype=float[8]):
plot(X,Y, color=COLOR(HUE, C), thickness=6);

The best method is probably using dsolve + fsolve, because fsolve + int(..., numeric) is too slow.

restart:
Digits:=15:
#eqList := ...
FS:=proc(a,f,t) 
    local u, ds, F;
    ds:=dsolve({diff(u(t),t) = f, u(-4)=0}, numeric, output=operator);
    F:=eval(u,ds);
    fsolve(F-a)
end proc: 
G:=e -> FS(op(1,e), op([2,1],e), op([2,2,1], e)):
map(G, eqList);

 [-4.00000000000000, -3.13995448471848, -2.11370056919619, 
  -1.52883183230736, -1.04840878317075, -0.706638188721254,
   -0.479577206139816, -0.315169354360080, -0.182432574734532, 
   -0.0625487819911634, 0.0891458203730475, 0.368690494498621, 
   1.50469145760957, 2.11552551065259, 2.46119669196488, 
   2.72088546517742, 2.96213514705916, 3.21307491249018, 
   3.59812548343646, 5.16349647450855, 6.56002211630892, 
   6.81413891022544, 6.66825093298453, 6.40577992538593, 
   6.19352683224088, 6.02201310955201, 5.95433562651451, 
   6.06761314123142, 6.32677648880848]

nops(eqList)=nops(%);
                            29 = 29

restart
L := [3*a0 + 2*a1, -a1/3 + a0/2 + 12*a3]:
x:=a3:

solve(select(has, L, x), x);

ex := x^4 + 3*(1+x^2)*f(x) + (x*x^2+x+1)*((D)(f)(x))^2 + (3*x+3)*(D@D)(f)(x):
eval( ex, f = (x -> (a1+a2*x+c*x^2+O(x^3)) ) ): 
series(%, x);
 #                   a2^2  + 3 a1 + 6 c + O(x)

Next time please post text, not pictures.

You have a polynomial function f in 2 variables and total degree 2, so
f(x,y) = a*x^2+b*y^2+c*x*y+p*x+q*y+r

You want that both partial functions f(., y0),  f(x0, .)   be concave [strictly].
This happens if and only if  a<0 and b<0.

Please note that f is concave (in both variables) iff a<0, b<0 and 4*a*b - c^2 > 0.

You post so many almost identical questions.
Here, the answer of your question why simplify does not work on this example? - MaplePrimes
applies ad litteram (but for x>0 now).

You have syntax erors. Use:

s := Maximize(TRC(tau1, lambda), C3 union C5, tau1 = 0 .. 1, lambda = 0 .. 1, assume = nonnegative);

What polynomials are congruent to u(x) modulo  p(x).

Answer: the equivalent class   < rem(u(x), p(x), x) > 
i.e.  rem(u(x), p(x), x) + A(x)*p(x),  where A(x) is an arbitrary polynomial,

Your example: 

rem(x^2+x+1, x+1, x) + A(x)*(x+1);  # for A(x)=x ==>  x^2+x+1

                               1 + A(x)*(x + 1)

In the last loop:
-  Insert a space between from and 1
and
-  replace sum with add.

Li is in the sense of Cauchy principal value, so, you must use:

int(1/ln(x),x=0..X, CPV);
#                           Ei(ln(X))

convert(%, Li);
#                             Li(X)

Just use := for assignments, instead of =.
So, 

n1 := 423; x := 16;   etc

Use

map(eval~, M);