Markiyan told you about the decision in the complex domain.
In the real domain, this equation has only two solutions  {a=0, b=0}  and  {a=1, b=1} . This is evident from the plot of your system:

sys := [b-a*sqrt(1+a^2+b^2)-a^2*(a*b-sqrt(1+a^2+b^2)) = 0, a-b*sqrt(1+a^2+b^2)-b^2*(a*b-sqrt(1+a^2+b^2)) = 0]:

plots[implicitplot](sys, a = -5 .. 5, b = -5 .. 5, color = [red, blue], numpoints = 50000);

For a rigorous solution try the command  RealDomain[solve]  or solve by hand. In the manual solution it is useful  firstly to factor the equations:

sys := [b-a*sqrt(1+a^2+b^2)-a^2*(a*b-sqrt(1+a^2+b^2)) = 0, a-b*sqrt(1+a^2+b^2)-b^2*(a*b-sqrt(1+a^2+b^2)) = 0]:

factor(sys);

Look this application  http://www.maplesoft.com/applications/view.aspx?SID=146470

restart;

for i do

a:=ithprime(i):

if a>300 then break fi:

L[i]:=a:

end do:

L:=[seq(L[k], k=1..i-1)];

L := [2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37, 41, 43, 47, 53, 59, 61, 67, 71, 73, 79, 83, 89, 97, 101, 103, 107, 109, 113, 127, 131, 137, 139, 149, 151, 157, 163, 167, 173, 179, 181, 191, 193, 197, 199, 211, 223, 227, 229, 233, 239, 241, 251, 257, 263, 269, 271, 277, 281, 283, 293]

Output is  more compact and can be used in the future, for example:

nops(L),  L[20];
                                   62, 71

Thanks Markiyan for your latest comment! This  hinted to me the idea of ​​solution.                                                    

First, we add the two equations. The resulting equation is a consequence of the original system and contains all of its roots. Then we estimate the left and right sides:

sys := [sqrt(sin(x)^2+1/sin(x)^2)+sqrt(cos(y)^2+1/cos(y)^2) = sqrt(20*y/(x+y)), sqrt(sin(y)^2+1/sin(y)^2)+sqrt(cos(x)^2+1/cos(x)^2) = sqrt(20*x/(x+y))]:

A := lhs(sys[1])+lhs(sys[2]);   B := rhs(sys[1])+rhs(sys[2]);

M := minimize(2*(eval(op(1, A)+op(2, A), [sin(x)^2 = p, 1/sin(x)^2 = 1/p, cos(x)^2 = 1-p, 1/cos(x)^2 = 1/(1-p)])), p = 0 .. 1, location);

N := maximize(eval(B, y = u*x), u = 0 .. infinity, location);

We see that the minimum value of the left-hand side is equal to the maximum value of the right side. The maximum value of the right-hand side is attained for any  x=y  (of cause  x<>0 and y<>0). The minimum value of the left-hand side is attained for any  x  and  y  such that  sin(x)^2 = 1/2  and  cos(x)^2 = 1/2 .

We find:

allvalues(solve({x = y, cos(x)^2 = 1/2, sin(x)^2 = 1/2}, AllSolutions));

The resulting solutions can be written more compactly:

x = Pi/4 + Pi*n/2,  y = Pi/4 + Pi*n/2,  n is integer

Since  the equation  A=B  is  a consequence of the original system, the solutions should be checked:

simplify(eval(sys, {x = (1/4)*Pi+(1/2)*n*Pi, y = (1/4)*Pi+(1/2)*n*Pi, 1/cos(x)^2 = 2, cos(x)^2 = 1/2, 1/cos(y)^2 = 2, cos(y)^2 = 1/2, 1/sin(x)^2 = 2, sin(x)^2 = 1/2, 1/sin(y)^2 = 2, sin(y)^2 = 1/2})) assuming n::integer;

                      [sqrt(5)*sqrt(2) = sqrt(5)*sqrt(2), sqrt(5)*sqrt(2) = sqrt(5)*sqrt(2)]

We prove that if the number  n  contains a lot of numbers (more than 6) it will not double.

At first, the idea of ​​proof by example. Let  n  has 7 digits. Then inequality  1000000 <= n <= 9999999  holds. We have 2000000<=2* n<=19999998 . Then  2*n[in base 7] <= 7^7+add(6*7^k, k=0..6) = 1647085 < 2000000

The general case reduces to the proof the inequality 

7^m+sum(6*7^k, k=0..m-1)<2*10^(m-1);   # m>=7

                        2*7^m - 1 < 2*10^(m-1)

Is sufficient to prove the stronger inequality  2*7^m < 2*10^(m-1)  for  m>=7  . It is equivalent to  10<(10/7)^m

isolve(10<(10/7)^m);
about(_NN1);

Checking for  n<1000000

t:=time():

N:=0:

for a0 from 0 to 1 do

for a1 from 0 to 6 do

for a2 from 0 to 6 do

for a3 from 0 to 6 do

for a4 from 0 to 6 do

for a5 from 0 to 6 do

for a6 from 0 by 2 to 6 do

a:=add(a||i*7^(6-i), i=0..6): b:=add(a||i*10^(6-i), i=0..6):

if 2*a=b then N:=N+1: L[N]:=a fi:

od: od: od: od: od: od: od:

[seq(L[i], i=1..N)];

time()-t;

                               [0, 51, 102, 105, 153, 156, 207, 210, 258, 261, 312, 315]
                                                                       2.188

Carl. your code is compact and elegant, but it works too slowly. Can you explain why?

Instead of  f(2)  write  eval(f(2), a = 2)  or change the procedure:

f := proc (b)

if b <= 0 then 0

elif b <= evalf(e*sin((1/2)*B)) then eval(eq1, a = b)

elif b <= evalf(2*e*sin((1/2)*B)) then eval(eq2, a = b)

else 0 end if

end proc;

Example:

f(2);

      1.358226754

restart;

g := proc(i)

if i = 1 then a else 0 fi

end proc:

h :=f->sum('g(i)', i=1 .. f):

h(3);

                 a

The use of the package  DirectSearch  is not proof in the mathematical sense.

The problem reduces to the optimization of a function of two variables because we can assume that  a=x, b=y, c=1  with restrictions

x^2<=y^2+1, y^2<=x^2+1, 1<=x^2+y^2

I took nonstrict inequalities, as maximum and minimum may be achieved within the domain or on its boundary.

Plot of the domain

plots[implicitplot]([x^2+y^2=1, y^2-x^2=1, x^2-y^2=1], x=0..4, y=0..4, color=black, thickness=2);

The domain of the function  is restricted by 3 smooth lines.

Find the function:

Expr:=proc(a,b,c)

local p, S, R;

p:=(a+b+c)/2;

S:=sqrt(p*(p-a)*(p-b)*(p-c));

R:=a*b*c/4/S;

simplify(R*p/(2*a*R+b*c));

end proc:

f:=unapply(Expr(x,y,1), x,y); 

Find the critical points of  f  within the domain:

solve({diff(f(x,y),x),diff(f(x,y),y)});

simplify(eval(f(x,y), %)); 

Find maximum and minimum of  f  on the boundaries of the domain: 

A:=subs([x=cos(t), y=sin(t)], f(x,y)):

simplify([maximize(A, t=0..Pi/2, location)]); evalf(%);

simplify([minimize(A, t=0..Pi/2, location)]); evalf(%);

B:=subs([x=cosh(t), y=sinh(t)], f(x,y)):

simplify([maximize(B, t=0..infinity,location)]); evalf(%);

simplify([minimize(B, t=0..infinity, location)]); evalf(%);

C:=subs([y=cosh(t), x=sinh(t)], f(x,y)):

simplify([maximize(C, t=0..infinity,location)]); evalf(%);

simplify([minimize(C, t=0..infinity, location)]); evalf(%);

The inequality  2/5 <= R*p/(2*a*R+b*c) < 1/2  is prooved. The lower limit  2/5  is reached for the triangle with sides  6/5*C, C, C, where C is arbitrary positive. The upper limit  1/2  is not achieved for any acute triangle.

Tuples:=proc(n, b)

local L, It;

L:=[seq([k], k=0..b)];

if n=1 then return L fi;

It:=proc(M)

[seq(seq([k, op(M[i])], k=0..b), i=1..nops(M))];

end proc;

(It@@(n-1))(L);

end proc;

Example:

Tuples(4, 2);

[[0, 0, 0, 0], [1, 0, 0, 0], [2, 0, 0, 0], [0, 1, 0, 0], [1, 1, 0, 0], [2, 1, 0, 0], [0, 2, 0, 0], [1, 2, 0, 0], [2, 2, 0, 0], [0, 0, 1, 0], [1, 0, 1, 0], [2, 0, 1, 0], [0, 1, 1, 0], [1, 1, 1, 0], [2, 1, 1, 0], [0, 2, 1, 0], [1, 2, 1, 0], [2, 2, 1, 0], [0, 0, 2, 0], [1, 0, 2, 0], [2, 0, 2, 0], [0, 1, 2, 0], [1, 1, 2, 0], [2, 1, 2, 0], [0, 2, 2, 0], [1, 2, 2, 0], [2, 2, 2, 0], [0, 0, 0, 1], [1, 0, 0, 1], [2, 0, 0, 1], [0, 1, 0, 1], [1, 1, 0, 1], [2, 1, 0, 1], [0, 2, 0, 1], [1, 2, 0, 1], [2, 2, 0, 1], [0, 0, 1, 1], [1, 0, 1, 1], [2, 0, 1, 1], [0, 1, 1, 1], [1, 1, 1, 1], [2, 1, 1, 1], [0, 2, 1, 1], [1, 2, 1, 1], [2, 2, 1, 1], [0, 0, 2, 1], [1, 0, 2, 1], [2, 0, 2, 1], [0, 1, 2, 1], [1, 1, 2, 1], [2, 1, 2, 1], [0, 2, 2, 1], [1, 2, 2, 1], [2, 2, 2, 1], [0, 0, 0, 2], [1, 0, 0, 2], [2, 0, 0, 2], [0, 1, 0, 2], [1, 1, 0, 2], [2, 1, 0, 2], [0, 2, 0, 2], [1, 2, 0, 2], [2, 2, 0, 2], [0, 0, 1, 2], [1, 0, 1, 2], [2, 0, 1, 2], [0, 1, 1, 2], [1, 1, 1, 2], [2, 1, 1, 2], [0, 2, 1, 2], [1, 2, 1, 2], [2, 2, 1, 2], [0, 0, 2, 2], [1, 0, 2, 2], [2, 0, 2, 2], [0, 1, 2, 2], [1, 1, 2, 2], [2, 1, 2, 2], [0, 2, 2, 2], [1, 2, 2, 2], [2, 2, 2, 2]]

Of cause, 25^2+24^2 . You can see it:

with(plottools):

A:=curve([[0,0],[200,0],[200,200],[0,200],[0,0]], thickness=2, color=black):

B:=seq(seq(disk([4+8*i,4+8*j], 1.5, color=green), j=0..24), i=0..24):

C:=seq(seq(disk([8+8*i,8+8*j], 1.5, color=yellow), j=0..23), i=0..23):

plots[display](A, B, C, axes=none);  # All in one

plots[display](A, B, axes=none);  # Only green trees

plots[display](A, C, axes=none);  # Only yellow trees

Two squares: the first square of green trees (with a side of 25 trees), and the second one of yellow trees (with the side of 24  tree):

I usually work in the classic interface. The text of the code can without any problems be copied and pasted into a text editor of mapleprimes. But if you do so, then the front of each line of code check mark appears. So at first I copy the code into Word, and then from Word into the text editor.

If I work in the standard interface, to copy the code I first select it, and then by context menu convert to 1-D Math, then cope without any problems .

for C from 2 to 10 do

s[C] := lhs(op(allvalues(solve({K > 0, K*(K-1) > 6*C-2})))):

end do:

L := [seq(floor(s[i]+1), i = 2 .. 10)];

                                      L := [4, 5, 6, 6, 7, 7, 8, 8, 9]

If you make a change  x=sqrt(lambda) , lambda>=0  it can be clearly seen from the graphs

plot([tan(x), x], x=-Pi..10*Pi, -5..35);

that in each range  x = (1/2)*Pi+Pi*k .. (1/2)*Pi+Pi*(k+1), k>=-1 there is a single root.

Finding  the first 10 roots:

seq(fsolve(tan(x) = x, x = (1/2)*Pi+Pi*k .. (1/2)*Pi+Pi*(k+1))^2, k = -1 .. 8);

0., 20.19072856, 59.67951595, 118.8998692, 197.8578111, 296.5544121, 414.9899843, 553.1646459, 711.0784498, 888.7314224

PS. This is interesting: my list does not coincide with Carl's one.

An interesting problem!

It can be solved in different ways. In my view, the most short way is to use a double integral with the change of variables. The change  u=y/x, v=(a-x)/y  maps the original region on the rectangle.

restart;

solve({u=y/x, v=(a-x)/y}, {x,y}):

assign(%):

int(Student[MultivariateCalculus][Jacobian]([x,y], [u,v], output = determinant), [u=1/2..1, v=1..3]) assuming a>0; 

                                                            7/120*a^2

A := plots[implicitplot](max(2-r, r-5, 3*Pi*(1/4)-theta, theta-5*Pi*(1/4)) = 0, r = 0 .. 6, theta = 0 .. 2*Pi, coords = polar, axiscoordinates = polar, gridrefine = 3):

B := plottools[polygon]([[2*cos(3*Pi*(1/4)), 2*sin(3*Pi*(1/4))], [5*cos(3*Pi*(1/4)), 5*sin(3*Pi*(1/4))], seq([5*cos(3*Pi*(1/4)+(1/200)*Pi*i), 5*sin(3*Pi*(1/4)+(1/200)*Pi*i)], i = 1 .. 100), [2*cos(5*Pi*(1/4)), 2*sin(5*Pi*(1/4))], seq([2*cos(5*Pi*(1/4)-(1/100)*Pi*i), 2*sin(5*Pi*(1/4)-(1/100)*Pi*i)], i = 1 .. 49)], color = green):

plots[display](A, B);

To automate the plotting of complicated plane figures, and to calculate their areas and perimeters, you can see my work

http://www.maplesoft.com/applications/view.aspx?SID=146470