@Kitonum Thank for your code. I tried
 

restart:
k := int(cos(3*x)^21, x):
algsubs(cos(3*x)^2 = 1-sin(3*x)^2, %):
a := factor(%):
b := expand(algsubs(sin(3*x) = t, a)):
algsubs(t = sin(3*x), %)

@acer Thank you very much.

@Kitonum How about tan(x)^7?
 

restart; 
`assuming`([int(tan(x)^n, x)], [n::posint]);
 eval(%, n = 9); 
expand(simplify(value(%))); 
algsubs(tan(x)^2+1 = 1/cos(x)^2, %)

I could not get the form

tan(x)^8/8 - tan(x)^6/6 + tan(x)^4/4 - tan(x)^2/2 - ln(cos(x))

@Kitonum When n is an odd number, n = 7, I like output in the form

restart; 
`assuming`([int(sin(x)^n, x)], [n::posint]):
eval(%, n = 7); 
expand(simplify(value(%))); int(sin(x)^7, x);

With cot(x)^20, I got the same form

@Kitonum I can't solve the first problem by myselt.

I considered tetrahedrons OABC with O(0,0,0), A(a, 0, 0), B(0, b, 0) and C(0, 0, c), where 0 < a < b < c<=100. 

@vv Note the second condition "The coordinates centre of sphere are also triples of integers".

@Kitonum Thank you very much.

@Kitonum 
How can I put like this?

Vertices([[[2, 5], [4, 3]], [[2, 5], [4, 9]], [[2, 6], [4, 8]], [[2, 7], [4, 
   3]], [[2, 7], [4, 5]], [[2, 7], [4, 9]]], -20 .. 20)

@Kitonum If I have a list[[[2, 5], [4, 3]], [[2, 5], [4, 9]], [[2, 6], [4, 8]], [[2, 7], [4, 
   3]], [[2, 7], [4, 5]], [[2, 7], [4, 9]]], where with the first element, [2, 5] is a centroid and  [4, 3] is orthocenter. How can I use this code?

@acer Thank you very much.

@acer How to select the number a and b so that the system of equations (x+1)^2+(y+3)^2 = 125  and (x-a)^2+(y-b)^2 = 225 have two integral solutions.

@Kitonum How can I use this code in geometry 2 D?

@vv Thank you very much.

@Carl Love Is there a procedure for every equation?

@Kitonum Thank you very much