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

x1^4 + x2^4 + 0.4*sin(7*x1) + 0.3*sin(4*Pi*x2) -1 = 0;
The task itself is not needed by anyone, it's just for rest and distraction from work.
(For example, this is not difficult to do with Draghilev's method.)

When I work with the fsolve and the number of variables is greater than 9, I have to use the op very carefully, because the order of the variables changes. The same thing happens when working with dsolve(numeric), which is especially inconvenient. For example, the RootFinding [Isolate]  works fine.
I ran into this a very long time ago and I have examples with 15 variables. It is completely unclear why such was done.
How to work with it or how to avoid it?
Example of the same solution using  fsolve and RootFinding [Isolate]

restart:
 CD1 := -.5; CD2 := 2; CD3 := 1.; 
 g1 := 5.; g2 := 3.; g3 := 2.; 
 cg1 := 0.; cg2 := 3.7; cg3 := 1;
 cd1 := 7.5; cd2 := 0; cd3 := 0; 
 L1 := .72; L2 := 8.6; L3 := 5.5; L4 := 1.25; L5 := 3.102; L6 := 5.1; L7 := 7.; 
f1 := (CD1-x4)^2+(CD2-x5)^2+(CD3-x6)^2-L1^2; 
f2 := x1-5; 
f3 := x5-1.45; 
f4 := (g1-x1)^2+(g2-x2)^2+(g3-x3)^2-L3^2; 
f5 := (x7-x1)^2+(x8-x2)^2+(x9-x3)^2-L2^2; 
f6 := (cg1-x7)^2+(cg2-x8)^2-(cg3-x9)^2-L4^2; 
f7 := x7+x8-1.2*cg2; 
f8 := (x4-x7)^2+(x5-x8)^2+(x6-x9)^2-L5^2; 
f9 := (cd1-x10)^2+(cd2-x11)^2+(cd3-x12)^2-L6^2; 
f10 := (x1-x10)^2+(x2-x11)^2+(x3-x12)^2-L7^2;
f11 := x11-.1*x12; 
f12 := x1-x2-x3-x4+x5+x6+x7+x8+x9+x10+x11-x12; 
fsolve({f1, f2, f3, f4, f5, f6, f7, f8, f9, f10, f11, f12}); 
RootFinding[Isolate]([f1, f2, f3, f4, f5, f6, f7, f8, f9, f10, f11, f12], [x1, x2, x3, x4, x5, x6, x7, x8, x9, x10, x11, x12])[4];
  {x1 = 5.000000000, x10 = 6.265979021, x11 = 0.4923895285, 
  x12 = 4.923895285, x2 = 7.363909611, x3 = 5.347580157, 
    x4 = -0.4802494112, x5 = 1.450000000, x6 = 0.5357695462,
    x7 = 2.014012007, x8 = 2.425987993, x9 = -1.029002454}
[x1 = 5., x2 = 7.363909611, x3 = 5.347580157, x4 = -0.4802494112, 
  x5 = 1.450000000, x6 = 0.5357695462, x7 = 2.014012007,
  x8 = 2.425987993, x9 = -1.029002454, x10 = 6.265979021,
  x11 = 0.4923895285, x12 = 4.923895285]