restart; with(geom3d): with(linalg):

A:=<0,0,0>: # A on l

B:= <2,-5,6>: # B on l

M :=<1,-1,3>:

P:=a*(x-A[1])+b*(y-A[2])+c*(z-A[3]):

Sol:=[solve({subs(x=B[1],y=B[2],z=B[3],P)=0,a^2+b^2+c^2=1,

abs(subs(x=M[1],y=M[2],z=M[3],P))=1},{a,b,c})]: # Equation of planes passing l and has a distance equal 1 from M

Equation(plane(alpha,subs(Sol[1], P)=0,[x,y,z])):

Equation(plane(beta,subs(Sol[2], P)=0,[x,y,z])):

n1:=NormalVector(alpha):

n2:=NormalVector(beta):

plane(Q,x + y + 5*z = 0,[x,y,z]): #Equation Q passing M and parallel to P

n:=NormalVector(Q):

s:= crossprod(n1,n):

r:=crossprod(n2,n):

convert(M,list):

point('M',op(%)):

Equation(line(Delta_1,[M,s],t));

Equation(line(Delta1_2,[M,r],t));

> restart; with(geom3d): with(LinearAlgebra):

A:=<0,0,0>:

B:= <2,-5,6>:

M :=<1,-1,3>:

P:=a*(x-A[1])+b*(y-A[2])+c*(z-A[3]):

Sol:=[solve({subs(x=B[1],y=B[2],z=B[3],P)=0,a^2+b^2+c^2=1,

abs(subs(x=M[1],y=M[2],z=M[3],P))=1},{a,b,c})]:

Equation(plane(alpha,subs(Sol[1], P)=0,[x,y,z])):

Equation(plane(beta,subs(Sol[2], P)=0,[x,y,z])):

n1:=convert(NormalVector(alpha), Vector):

n2:=convert(NormalVector(beta),Vector): plane(Q,x + y + 5*z = 0,[x,y,z]): n:=convert(NormalVector(Q),Vector):

s:= CrossProduct(n1,n): r:=CrossProduct(n2,n):

N:=<x,y,z>: [seq(N[i]=(M + s*t)[i],i=1..3)];

[seq(N[i]=(M + r*t)[i],i=1..3)];

restart;

with(geom3d):

line(d,[2*t + 3, t-2, -t-1],t):

plane(P,x + y + z + 2 = 0,[x,y,z]):

coordinates(intersection(M,d,P)):

with(LinearAlgebra):

H:=<x,y,z>:

a:=convert(ParallelVector(d),Vector):

n:=convert(NormalVector(P),Vector):

v:=CrossProduct(a,n):

m:=convert(coordinates(M),Vector):

Sol:=solve([x + y + z + 2 = 0, DotProduct(m - H, v,conjugate=false)=0, Norm(m - H,2)^2 = 42]):

H1:=[seq(rhs(Sol[1,i]),i=1..3)]:

H2:=[seq(rhs(Sol[2,i]),i=1..3)]:

b:=convert(H1,Vector):

c:=convert(H2,Vector):

 A:=<x,y,z>:

[seq(A[i]=(b + v*t)[i],i=1..3)];

[seq(A[i]=(c + v*t)[i],i=1..3)];

solve([x + y + z + 2 = 0, DotProduct(m - H, v,conjugate=false)=0, Norm(m - H,2)^2 = 42]);

Please help me. Thank you very much.

 restart;

with(geom3d):

line(d,[2*t + 3, t-2, -t-1],t):

plane(P,x + y + z + 2 = 0,[x,y,z]):

coordinates(intersection(M,d,P)):

with(LinearAlgebra):

H:=<x,y,z>:

a:=convert(ParallelVector(d),Vector)

n:=convert(NormalVector(P),Vector):

v:=CrossProduct(a,n):

m:=convert(coordinates(M),Vector);

solve([x + y + z + 2 = 0, DotProduct(m - H, v,conjugate=false)=0, Norm(m - H,2)^2 = 42],[x,y,z]);

Please consider for me at the end line. Thank you very much.

This is second my code.

restart;with(geom3d):

plane(P,3*x+12*y-3*z-5=0,[x,y,z]):

plane(Q,3*x - 4*y  + 9*z + 7=0,[x,y,z]):

line(d1, [2*t -5, -4*t+3, 3*t -1],t):

line(d2, [-2*m + 3, 3*m - 1, 4*m + 2],m):

with(LinearAlgebra):

a:=convert(NormalVector(P),Vector):

b:=convert(NormalVector(Q),Vector):

v:=CrossProduct(a,b):

A:=<2*t -5, -4*t+3, 3*t -1>:

B:=<-2*m+3, 3*m - 1, 4*m + 2>:

u:=A - B:

w:=LinearAlgebra[CrossProduct](v,u):

solve([seq(w[i]=0,i=1..3)]): assign(%):

'A'=A: 'B'=B:

M:=<x,y,z>: [seq(M[i]=(A + v*n)[i],i=1..3)];

Please comment for me. Thank you very much.

Another way.

with(geom3d):
plane(P,3*x+12*y-3*z-5=0,[x,y,z]):
plane(Q,3*x - 4*y  + 9*z + 7=0,[x,y,z]):
line(d1, [2*t -5, -4*t+3, 3*t -1],t):
line(d2, [-2*m + 3, 3*m - 1, 4*m + 2],m):
intersection(Delta,P,Q):
delta:=ParallelVector(Delta):
a:=ParallelVector(d1):
with(linalg):
n1:= crossprod(delta,a):
plane(P1,[point(A,-5,3,-1),n1],[x,y,z]):
coordinates(intersection(M,d2,P1)):
line(d,[M,delta],t):
Equation(d);
AreParallel(d,d1);
AreParallel(d,d2);
                  
                            

I am sory about that. 

Let a = (2,1,-1) be direction of (d);

n = (1,1,1) be normal vector of (P);

v  be direction of Delta. Then v = CrossProduct(a, n). 

The statemnet 2) need to repair  Vector(MH) perpendicular to vector v.

Note that M(1,-3,0) is intersection of point of (d) and (P).

restart:

with(geom3d):

line(d1,[1+t, -2-t,t],t):

line(d2,[1 -m, 3-2*m,m - 1],m):

a:=convert(ParallelVector(d1), Vector):

b:=convert(ParallelVector(d2), Vector):

with(LinearAlgebra):

A:=<1+t, -2-t,t>:

B:=<1 -m, 3-2*m,m - 1>:

sys:=solve([DotProduct(A-B, a,conjugate=false) = 0, DotProduct(A-B, b,conjugate=false) = 0],[t,m]): assign(%):

'A'=A: 'B'=B: w:=A-B:

M:=<x,y,z>: [seq(M[i]=(A + w*n)[i],i=1..3)];

>AreSkewLines(d1,d2);

This is my idea: Let (P) be finding plane.  Put n = (a, b, c) is a normal vector of (P). Then, the equation of (P) has the form

a*(x - 1) + b*(y -3) + c(z + 6) = 0. We find a, b, c from following conditions:

1) a^2 + b^2 +c^2 = 1;

2) vector n is perpendicular to vector AB;

3) The distance form C(1; -1; 7) (C is center of the given sphere) to the plane (P) equal to 3.

But i assume(a^2 + b^2 + c^2 > 0) and use command 

plane(P,[A, n],[x, y, z]) can not write the Equation of (P). Plese help me. Thank you very much.

I used your code. 

restart:
with(LinearAlgebra):A:=<-1,3,6>: B:=<2,2,-0>: N:=<x,y,z>: o:=<1,-1,7>:
{Norm(o-N,2)=3, DotProduct(N - o, N-A) = 0, DotProduct(N - o, N-B) = 0}:
solve({Norm(o-N,2)=3, DotProduct(N - o, N-A) = 0,DotProduct(N - o, N-B) = 0}): assign(%): N:
'N'=[seq(N[i],i=1..3)];

and i got Warning, solutions may have been lost.

Please check my code. Thank you.

with(LinearAlgebra):

A:=<0,1,2>:  M:=<2*t,1 + t,-1 - t>: N:=<1+m, -1 -2*m,2+m>:

u:=A - M: v:=A-N:
w:=CrossProduct(u,v);

sys:=solve([t*m-(3+t)*(2+2*m) = 0, (3+t)*(-1-m) -2*t*m=0, 2*t*(2+2*m) +t*(-1-m)=0],[t,m]);

subs(sys[1],M);

subs(sys[1],N);

In this way, i don't know how to solve w = 0 in order to find m and t automatically. Please help me. Thank you very much.