@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 

@Kitonum I can't get the result with this code.
 

restart;
k := 0;
for a from -5 to 15 do
for b from -5 to 15 do
for c from -5 to 15 do
if igcd(b, c) = 1 and nops({a, b, c}) = 3 then X := [solve(x = a+sqrt(b*x+c))];
if type(X[1], integer) and type(X[2], integer) then k := k+1; L[k] := [a, b, c, X[]] end if end if end do end do end do;
L := convert(L, list);
k;
L;

@Kitonum Thank you very much.

@Kitonum where I add the condition nops({a,b,c,d,e,f})= 6? And some results are not correct. For example
solve(x = (3*x^2+6*x+2)/(x^2+5*x+10), x)

Mathematica out put 37 solutions.