@Markiyan Hirnyk  Your code does not work for an equilateral triangle, because in this case, O2 = O3. In general, the code is very cumbersome , because to find the point   O1, O2, O3  and  etc. you actually use the same calculations that can be replaced by a single procedure.

@Markiyan Hirnyk  In your code you have problems with  xC := 1/2:  yC := (1/2)*sqrt(3):

@Markiyan Hirnyk   It is not interesting. To prove just replace  a= ln(4)+sin(1)

@Markiyan Hirnyk  We must use the exact math! Replace  1.1  with  11/10 .

There is a more interesting case, when the code is not working. Try to find it!

@Markiyan Hirnyk  Thank you for the detailed checking my code! Of course, if we consider only the proof of the original assertion, the proof itself conducted to the general case (for 6 points). If in the particular case  a few points are merged into one, then by continuity it follows that these points also lie on a circle.

Here is the code that works for all cases and even for degenerate triangle and  returns also the center of the circle and its radius.

restart;

Cc:=proc(A,B,C)

local x1, y1, x2, y2, x3, y3, x, y;

x1,y1:=op(A);  x2,y2:=op(B);  x3,y3:=op(C);

if x2*y3-x2*y1-x1*y3-x3*y2+x3*y1+x1*y2=0 then error "No triangle. Three points are collinear" fi;

solve({(x2-x1)*(x-(x1+x2)/2)+(y2-y1)*(y-(y1+y2)/2)=0, (x2-x3)*(x-(x2+x3)/2)+(y2-y3)*(y-(y2+y3)/2)=0},{x,y});

assign(%);

[simplify([x,y]), simplify(sqrt((x-x1)^2+(y-y1)^2))];

end proc:

The proof:

A, B, C:=[x1,y1], [x2,y2], [x3,y3]:

A1, B1, C1, M:=(B+C)/2, (A+C)/2, (A+B)/2, (A+B+C)/3:

P1:=Cc(A,M,B1)[1]: P2:=Cc(B1,M,C)[1]: P3:=Cc(C,M,A1)[1]:

P4:=Cc(A1,M,B)[1]: P5:=Cc(B,M,C1)[1]: P6:=Cc(C1,M,A)[1]:

P:={seq(P||i, i=1..6)}:

Cc1:=Cc(P[1],P[2],P[3]);

if is(eval((x-Cc1[1,1])^2+(y-Cc1[1,2])^2, {x=P4[1],y=P4[2]})=Cc1[2]^2) and

is(eval((x-Cc1[1,1])^2+(y-Cc1[1,2])^2, {x=P5[1],y=P5[2]})=Cc1[2]^2) and

is(eval((x-Cc1[1,1])^2+(y-Cc1[1,2])^2, {x=P6[1],y=P6[2]})=Cc1[2]^2) then true else false fi;

Why do you ask the same question?

http://www.mapleprimes.com/questions/201010-How-To-Computer-GSO-In-Maple#answer203907

@Carl Love   Of course, you're right! It is my carelessness.

@toandhsp   

eq1:=(x+y)*(x^2+y^2) = 5500;

eq2:=(x-y)*(x^2-y^2) = 352;

solve[RealDomain]({eq1,eq2});

@Markiyan Hirnyk  Maple does not compute this integral, since the exact form of the term  O(x - Pi)  as a function of  x  is unknown. But it is clear that this function is limited, and integral of limited function over a finite interval is limited. Therefore, it is obvious even without Maple that the sum of infinitely large and  limited functions is infinitely large function.

I will give this example to my students. Markiyan, please tell us about the source of this problem.

@Markiyan Hirnyk   These commands, that you listed, allow you to find the roots of the equation in the specified interval, but say nothing about the structure of the set of all solutions of this equation, for example, whether this set is finite or infinite, etc. In my post revealed the structure of this set of solutions that allows you to find any root  by its number. For example, it is easy to find a millionth root without calculating all the preceding roots.

@Axel Vogt  Thanks for the useful comments. Of course, equality  

x^t/y^t = (x/y)^t

is not always true. However, if  x>0  and  y>0  then it is true. 

Another interesting example:

1=sqrt((-1)*(-1))=sqrt(-1)*sqrt(-1)=I*I=I^2=-1

@Carl Love   Thank you for your detailed explanation!

@Carl Love  Thank you for your helpful comment! The effectiveness of your package (0.1 seconds) for solving such problems is striking. But you wrote that your method does not directly prove the uniqueness of the solutions found. And how much time require such proof?

@Markiyan Hirnyk  Well, and if  t = 10^8 .. 10^8+20 ?