The question is not at all from me, but, probably, one might say, from the authors of this publication.  interesting_system.pdf

Just for fun.
Find all real solutions to this 2x2 system of nonlinear equations in any given domain. 

f1:=x1-x1*sin(x1+5x2)-x2*cos(5x1-x2);
f2:=x2-x2*sin(5x1-3x2)+x1*cos(3x1+5x2);

Inscribed square problem

I decided to check on this curve 
 

 4*(x1-0.25)^4-x1^2*x2^2+(x2-0.25)^4-1.21=0;

I get a very rough solution, because the difference between the sides of the "square" begins already at 1-2 decimal places. More precisely, it doesn’t work, that is, we can say that I personally could not find confirmation of the hypothesis.
The coordinates of the vertices of the square:

                    -0.4823584672, -0.2770841741

                    0.9883885535, -0.3959790155

                    1.108267478, 1.086941264

                   -0.3459185869, 1.219514527

Side lengths:
                          1.475544911

                          1.487757882

                          1.460216690

                          1.502805215


 

Perhaps someone would like to try.

Another training example (number 2 and last) for finding all solutions to a system of equations:

f1 := x3^2-0.1*x1^4-0.05*x2^4+1;
f2 := x1^3+x2^3+0.05*x3^3-1; 
f3 := -2*cos(3*x1)+2*cos(3*x2)-2*cos(3*x3)+1;

In my version, there are 116 solutions.
Is it so?

This is a training example for finding all solutions to a system of equations. If you look at the graph, you can count 36 solutions, but I managed to find 20 relatively good approximations. No attempts to get more solutions, which are also visible as intersections of graphs, did not lead to success. Therefore, there is a suspicion that there are only 20 solutions.
Is it so?

 restart: with(plots):
 a:=8.:
 f1 := x1^4-1.999*x1^2*x2^2+x2^4-1; 
 f2 := tan(x1+x2)-x2*sin(x1);
 implicitplot([f1, f2], x1 = -a .. a, x2 = -a .. a, numpoints = 25000, scaling = constrained,  color = [red, blue], thickness = 1);

   1, (3.192246883291975), (-3.0395187374365404)
   2, (3.0952031367176476), (-3.2447717313041897)
   3, (0.5881900748267959), (-1.160066226905079)
   4, (-0.936866718243322), (-1.3700058362814254)
   5, (-2.555853694651265), (-2.7399958564861953)
  6, (-3.2556241416421168), (-3.3964651254113774)
   7, (-3.583319843955091), (-3.7077839724189228)
   8, (-5.364827188794712), (-5.401998918608201)
   9, (-5.398295356665546), (-5.360818510223991)
  10, (-3.769206506106412), (-3.6477855329362683)
  11,(-1.3978806247566642), (-0.9772190664843745)
   12, (-1.192159295544978),(0.6492335177657542)
   13, (-3.0867255059416623),(2.927375855548188)
    14, (-3.18519036357835), (3.329801919022179)
   15, (2.0108268901120754), (2.243492422396739)
   16, (3.1133812329649766), (3.261937603184373)
   17, (4.0265558604742715), (4.130826167761226)
    18, (4.124539552922121), (4.019977762680433)
    19, (3.172338340844501), (3.018365965761908)
   20, (2.1945695320368097), (1.9558412553082192)

A simple tutorial task on optimization from one of the forums. We need to find the minimum of the function 
g1 = (x1-x2)*(x2-x3)*(x3-x4)*(x4-x1) with the following constraint g2 = x1^2+x2^2+x3^2+x4^2-1. 
It reduces to solving polynomial equations.  
Minimum = - 0.25. It's  solution:
[x1 = .683012702, x2 = .183012702, x3 = -.183012702, x4 = -.683012702],
[x1 = .183012702, x2 = .683012702, x3 = -.683012702, x4 = -.183012702],
[x1 = -.183012702, x2 = .183012702, x3 = .683012702, x4 = -.683012702],
[x1 = -.683012702, x2 = .683012702, x3 = .183012702, x4 = -.183012702],
[[x1 = .683012702, x2 = -.683012702, x3 = -.183012702, x4 = .183012702],
[x1 = .183012702, x2 = -.183012702, x3 = -.683012702, x4 = .683012702],
[x1 = -.183012702, x2 = -.683012702, x3 = .683012702, x4 = .183012702],
[x1 = -.683012702, x2 = -.183012702, x3 = .183012702, x4 = .683012702]
when x5 =0 .25

But so far I have not found an easy way to solve it using Maple. First we need to be able to successfully use fsolve, and only then "polynomial" functions.

restart; with(RootFinding):
g1 := (x1-x2)*(x2-x3)*(x3-x4)*(x4-x1); 
g2 := x1^2+x2^2+x3^2+x4^2-1;
 g :=g1+x5*g2; 
f1 := diff(g, x1); 
f2 := diff(g, x2); 
f3 := diff(g, x3); 
f4 := diff(g, x4); 
f5 := diff(g, x5); 
#solve([f1, f2, f3, f4, f5], [x1, x2, x3, x4, x5]);
Isolate([f1, f2, f3, f4, f5], [x1, x2, x3, x4, x5]); 
S := fsolve([f1, f2, f3, f4, f5], {x1, x2, x3, x4, x5}, maxsols = 8); 
x5 := rhs(op(5, S));
Isolate([f1, f2, f3, f4], [x1, x2, x3, x4]); 
solve([f1, f2, f3, f4], [x1, x2, x3, x4])

The optimization package doesn't help much either.
The task itself is of no interest, it is interesting to look at its simplest solution in Maple.
(By the way, Draghilev's method works well, but, of course, this is not the easiest way).