Ax:=t->cos(5*t): Ay:=t->0: Bx:=t->-cos(5*t)/2: By:=t->3*cos(5*t)/2:

S:=seq(plots[arrow]([Ax(i*Pi/72),Ay(i*Pi/72)],color=blue),i=0..72):

T:=seq(plots[arrow]([Bx(i*Pi/72),By(i*Pi/72)],color=green),i=0..72):

U:=seq(plots[display]([S[i],T[i]]),i=1..73):

plots[display](U,view=[-2..2,-2..2],insequence=true);

You must assign  z  as a function, such as in an example

restart:

dsolve({diff(z(t),t) = 2*t-1, z(0) = 1});

z:= unapply(rhs(%), t);  

z(60);

a:= 1: b:= 1/sqrt(3):

for n to 2010 do

c:= simplify((a + b)/(1-a*b)):  a:=b:  b:=c:

od:

c;

Brute force computation shows that there are no solutions. Explanation of code: firstly founded a list of matrices  H  of 0 and 1, satisfying the first condition. The number of such matrices will be 3190. As an example displayed 4 such matrices. Then, all of these 3190 matrices tested for the second condition.

st:= time():

It:=proc(K)

[seq([0,op(K[i])],i=1..nops(K)),seq([1,op(K[i])],i=1..nops(K))];

end proc:

M:=(It@@5)([[ ]]):

H:=[]:

for L[1] in M do

for L[2] in M do

for L[3] in M do

for L[4] in M do

for L[5] in M do

if add(add(L[i][j],i=1..3),j=1..3)=5 and add(add(L[i][j],i=1..3),j=2..4)=5 and add(add(L[i][j],i=1..3),j=3..5)=5 and

   add(add(L[i][j],i=2..4),j=1..3)=5 and add(add(L[i][j],i=2..4),j=2..4)=5 and add(add(L[i][j],i=2..4),j=3..5)=5 and

   add(add(L[i][j],i=3..5),j=1..3)=5 and add(add(L[i][j],i=3..5),j=2..4)=5 and add(add(L[i][j],i=3..5),j=3..5)=5 then

   H:=[op(H), Matrix([L[1],L[2],L[3],L[4],L[5]])]: fi:

od: od: od: od: od:

'nops(H)'=nops(H);  seq(H[n],n=1001..1004);

P:=[ ]:

for k in H do

if add(add(k[i,j],i=1..2),j=1..4)=4 and add(add(k[i,j],i=1..2),j=2..5)=4 and add(add(k[i,j],i=2..3),j=1..4)=4 and

   add(add(k[i,j],i=2..3),j=2..5)=4 and add(add(k[i,j],i=3..4),j=1..4)=4 and add(add(k[i,j],i=3..4),j=2..5)=4 and

   add(add(k[i,j],i=4..5),j=1..4)=4 and add(add(k[i,j],i=4..5),j=2..5)=4 and add(add(k[i,j],j=1..2),i=1..4)=4 and

   add(add(k[i,j],j=1..2),i=1..4)=4 and add(add(k[i,j],j=1..2),i=1..4)=4 and add(add(k[i,j],j=1..2),i=1..4)=4 and

   add(add(k[i,j],j=1..2),i=1..4)=4 and add(add(k[i,j],j=1..2),i=1..4)=4 and add(add(k[i,j],j=1..2),i=1..4)=4 and

   add(add(k[i,j],j=1..2),i=1..4)=4

   then P:=[op(P), k]: fi:

od:

P;   

time() - st;

The evalf command after solve command for systems as well as fsolve gives only 1 solution.

Possible decision of this problem:

f:=x->x^3-3*x^2+2:

g:=x->k*(x+1)+3:

RealDomain[solve]([f(x) = g(x), diff(f(x),x) = diff(g(x),x)]);

evalf(allvalues([%]));

Got all the solutions in both symbolic and numerical form!

Your solution is right.

I suggest another variant to solve your problem without solve command by direct computation (firstly find the points of tangancy):

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

line(d1,[-5+2*t,-3*t+1,-13+2*t],t):  line(d2,[-7+3*t,-1-2*t,8],t):  

a:=ParallelVector(d1):  b:=ParallelVector(d2):  

c:=crossprod(a,b):

sphere(S, x^2+  y^2 + z^2 -10*x +2*y +26*z -113=0, [x,y,z]):

T:=coordinates(center(S)):  R:=radius(S): 

H:=[T + [seq(c[i]/norm(c,2)*R,i=1..3)], T - [seq(c[i]/norm(c,2)*R,i=1..3)]]:

for i to 2 do

sort(Equation(plane(P,c[1]*(x-H[i,1]) +c[2]*(y-H[i,2]) +c[3]*(z-H[i,3]) =0,[x,y,z])), {x,y,z}):

od;

I have not found a suitable way to sort, so I suggest a way to manually isolate complete squares.

The code:

restart:  with(geom3d):

a:=[-t, -t, t]: point(T, a): point(A, 1, -1, -1): point(B, 2, 1, 2): point(C, 1, 3, 1): plane(ABC, [A,B,C], [x, y, z]):

d:= distance(T, ABC):  R:= distance(T, A):

r:= R^2 - d^2;  simplify(r) assuming real;

sol:= minimize(r, location);  Sol:= op(sol[2])[1];

point(T, subs(Sol, a)):  R:= subs(Sol, R):

Eq:=lhs(Equation(sphere(S, [T,R], [x, y, z])));

(x+coeff(Eq,x)/2)^2+(y+coeff(Eq,y)/2)^2+(z+coeff(Eq,z)/2)^2=sqrt((coeff(Eq,x)^2+coeff(Eq,y)^2+coeff(Eq,z)^2)/4-tcoeff(Eq))^`2`;  # Canonical equation of S

You can use the parametric equations of the surface.

The example:

theta:=Pi/3:

plot3d( [u*cos(t)*cos(theta), u*cos(t)*sin(theta), u*sin(t)], u=0..1, t=0..2*Pi, color=red, axes=normal, scaling=constrained, style=surface, view=[-1..1, -1..1, -1.2..1.2] );

For example:

Sol:= 9.356887917*10^5 , { [{P = 330.1448618, W = 3.339008977}, 9.356887917*10^5 ] }:

P:=rhs(Sol[2,1,1,1]);  W:=rhs(Sol[2,1,1,2]);

P:=330.1448618

W:= 3.339008977

I would write your code so:

restart:  with(geom3d):

a:=[-t, -t, t]: point(T, a): point(A, 1, -1, -1): point(B, 2, 1, 2): point(C, 1, 3, 1):

plane(ABC, [A,B,C], [x, y, z]):

d:= distance(T, ABC):  R:= distance(T, A):

r:= R^2 - d^2;  simplify(r) assuming real;

sol:= minimize(r, location);  Sol:= op(sol[2])[1];

point(T, subs(Sol, a)):  R:= subs(Sol, R):

Equation(sphere(S, [T,R], [x, y, z]));

Eq:=Student[Precalculus][CompleteSquare](Equation(S));

op(1, Eq) - op(4, op(1, Eq)) = sqrt(-op(4, op(1, Eq)))^`2`;  # Canonical equation of S

Here you will find answers to your questions!

Sorry! I made a slip of the pen. Instead  2  must be  d .

After  v:=crossprod(w,a)  insert two lines into your code:

d:=igcd(v[1],v[2],v[3]):

v:=map(x->x/2,v):

The lower boundary your region is a common segment  [-1, 1] . So look for  maximum and minimum of the function of one variable  f(x,0)=x^2+x  on this segment by the commands  maximize  and  minimize . In this decision, the derivatives are not needed.

Computed plane passes through the line AB so that AB and its projection onto given plane  are perpendicular to the line of intersection of computed plane with given plane. Draw a picture and everything will be clear!

restart:

with(LinearAlgebra):

A:=<1, -2, 4>:  B:=<3, 5, -1>:  n1:=<1, 1, 1>:  M:=<x, y, z>:

n2:=CrossProduct(n1, B - A):  n:=CrossProduct(n2, B - A):  d:=igcd(seq(n[i],i=1..3)):

sort((n/d).(M - A)) = 0;   # Equation of the required plane