I am trying to find coodinates projection H(x,y, z)  of the point P :=<x3, y3, z3> on the line AB, where A:= <x1,y1,z1>, B:= <x2,y2,z2>. I tried

restart:
with(LinearAlgebra); 
A:= <x1,y1,z1>:
B:= <x2,y2,z2>:
P :=<x3, y3, z3>:
H:=<x,y,z>;
solve([DotProduct(B-A, H-P, conjugate = false) = 0, B-A = k*(H-A)], [x, y, z, k]);

I can not get the result. How can I get the coordinates of the point H?

 

Tetrahedron with vertices A(1,2,3), B(-2, 8, 9) , C(5, 0, 7) and D(3, 4, 2) has centre of sphere is inscribed I(2, 3, 4) and radius is r = 1. 
I tried
 

restart:
with(geom3d):
point(P1, 1, 2, 3):
point(P2, -2, 8, 9):
point(P3, 5, 0, 7):
point(P4, 3, 4, 2):
s1 := area(triangle(T1, [P2, P3, P4])):
s2 := area(triangle(T2, [P1, P3, P4])):
s3 := area(triangle(T3, [P1, P2, P4])):
s4 := area(triangle(T4, [P1, P2, P3])):
s := s1+s2+s3+s4:
x_I := simplify((xcoord(P1)*s1+xcoord(P2)*s2+xcoord(P3)*s3+xcoord(P4)*s4)/s):
y_I := simplify((ycoord(P1)*s1+ycoord(P2)*s2+ycoord(P3)*s3+ycoord(P4)*s4)/s):
z_I := simplify((zcoord(P1)*s1+zcoord(P2)*s2+zcoord(P3)*s3+zcoord(P4)*s4)/s):
plane(P, [P2, P3, P4], [x, y, z]):
point(T, x_I, y_I, z_I);
r := distance(T, P)
                           

I am trying to find all tetrahedrons that  sphere is inscribed in a  tetrahedron satisfy the following conditions:

1. Each vertex (x, y, z) is a triples of integers such that -20 <= x, y, z <= 20.
2. The coordinates centre of sphere are also triples of integers.

How can I find? 

I am trying write a code for this question at here 

https://math.stackexchange.com/questions/170319/how-many-triangles-with-integral-side-lengths-are-possible-provided-their-perim/170325#170325

but I got the wrong answer. This is my code.
 

restart; 
L := []; 
for a to 18 do
 for b to 18 do 
for c to 18 do 
if a < b+c and b < a+c and c < b+a and a > abs(b-c) and b > abs(a-c) and c > abs(a-b) and a+b+c = 36 then
 L := [op(L), [a, b, c]] end if end do end do end do;
 nops(L); 
L

I got 136 triangles. This is a wrong answer. How can I repair my code?

I am trying to find the number a and b (-20 < a < 20, -20 < b <20) so that two circles (x+1)^2+(y+3)^2 = 125 and (x-a)^2+(y-b)^2 = 225 cut at two points A and B and coordinates A and B are pairs of integers.
I tried
 

restart; L := []; 
for a from -20 to 20 do 
for b from -20 to 20 do 
for x from -20 to 20 do 
for y from -20 to 20 do 
if (x+1)^2+(y+3)^2 = 125 and (x-a)^2+(y-b)^2 = 225 and nops({a, b, x, y}) = 4 and x*y*a*b <> 0 then 
L := {op(L), {[a, b], [x, y]}} 
od: od: od: od:
nops(L); 
L

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. For example
 

restart; 
solve({(x-6)^2+(y+2)^2 = 225, (x+1)^2+(y+3)^2 = 125}, {x, y})

I want to find general formulas of the sequence satifies:
f(1) = 1, f(n+1) = 1/2*(f(n)+9/f(n))
I tried 
restart;
rsolve({f(1) = 1, f(n+1) = 1/2*(f(n)+9/f(n))}, {f})
I got the answer
3*coth((1/2)*arccoth(1/3)*2^n)

How to get the form (-3*2^(2^(-1+n))+3*(-1)^(1+2^(-1+n)))/(-2^(2^(-1+n))+1)?