You do some calculation, but the result of that calculation is not assigned to anything.
A similar example:

s:=2*x+3;
eval(s, x=1);
s;

Maple shows in detail all the steps only for solving linear equations.
Example:

restart;
Student:-Basics:-LinearSolveSteps((x+1)/a = 4/b+3*x/c, x);

There is a resource (not Maple) in which a very wide range of mathematical problems is solved step by step.

See  https://www.universalmathsolver.com/

The proof, that the trajectory of motion of point C  for any  L>0  when changing the parameter  a is a hyperbola ,  is easy to carry out using the  geometry package. We just denote the coordinates of point  C  by  x  and  y   and applying the condition  BC=CD  we get the equation of a hyperbola. For clarity, I also found the asymptotes of this hyperbola:

restart;
local D:
A:=[-L,0]: B:=[L,0]: C:=[x,y]: D:=[-x,y]:
Dist:=(X,Y)->sqrt((X[1]-Y[1])^2+(X[2]-Y[2])^2):
Eq:=(Dist(C,D)=Dist(C,B))^2;

# The proof

with(geometry):
conic(p, Eq, [x,y]) assuming L>0;
detail(p);
asymptotes(p);
y_acymp:=solve~(Equation~(asymptotes(p)), y);
y:=solve(Eq,y)[1];

# The procedure for plotting

P:=proc(X, L0)
local Curve, Asymptote, Trapezoid, T;
uses plots, plottools, geometry;
Curve:=plot(eval([y,-y],L=L0), x=-L0/3..15, color=red, thickness=3);
Asymptotes:=plot(eval([-sqrt(3)*x-(1/3)*sqrt(3)*L, sqrt(3)*x+(1/3)*sqrt(3)*L],L=L0), x=-L0/3..15, linestyle=3, color=black, thickness=0);
Trapezoid:=polygon(eval([A,B,C,D],[L=L0,x=X]),color="LightGreen");
T:=textplot([[eval(A,[x=X,L=L0])[],"A"],[eval(B,[x=X,L=L0])[],"B"],[eval(C,[x=X,L=L0])[],"C"],[eval(D,[x=X,L=L0])[],"D"]], align={above,left}, font=[TIMES,16]);
display(Curve, Asymptotes, Trapezoid, T, scaling=constrained, size=[400,800]);
end proc:

# Example:

a:=7;
P(a/2, 6);

Animation:

plots:-animate(P,['a'/2,6], 'a'=4..28, frames=90, size=[400,800]);

Not every tetrahedron has such a sphere. Such a sphere exists if and only if one of the following conditions is met:
a) the sums of the lengths of pairs of opposite edges are equal;
b) the circles inscribed in the faces touch in pairs;
c) the perpendiculars to the planes of the faces, raised from the centers of the circles inscribed in the faces, intersect at one point.

Consider a tetrahedron  ABCE  for which the edge lengths are known: AB=5, BC=6, AC=7, AE=5, BE=4, CE=6  (it is obvious that here the condition a) holds.) . We will assume that points A, B, C lie in the plane  xOy :  A(0,0,0), B(x,y,0), C(7,0,0) , E(u,v,w) . Next, we find the unknown coordinates of all the points, find the center  CC(x0,y0,z0)  and the radius R  of this sphere :

restart;
A:=[0,0,0]: B:=[x,y,0]: C:=[7,0,0]: E:=[u,v,w]:
Dist:=(X,Y)->sqrt((X[1]-Y[1])^2+(X[2]-Y[2])^2+(X[3]-Y[3])^2):
Dist1:=proc(X,Y,Z)
local XY, XZ, YZ, p, S;
XY:=Dist(X,Y); XZ:=Dist(X,Z); YZ:=Dist(Y,Z);
p:=(XY+XZ+YZ)/2; S:=sqrt(p*(p-XY)*(p-XZ)*(p-YZ));
2*S/Dist(Y,Z);
end proc:
solve({Dist(A,B)=5,Dist(B,C)=6,Dist(A,E)=5,Dist(B,E)=4,Dist(C,E)=6},explicit);
Sol:=select(s->`and`(seq(is(rhs(p)>=0),p=s)), [%])[];
assign(%):
CC:=[x0,y0,z0]:
solve({Dist1(CC,A,B)=R,Dist1(CC,B,C)=R,Dist1(CC,A,C)=R,Dist1(CC,A,E)=R,Dist1(CC,B,E)=R,Dist1(CC,C,E)=R}, explicit):
Sol1:=expand(%);
assign(Sol1):

with(plots): with(plottools):
AB:=line(A,B,color=red,thickness=2): BC:=line(B,C,color=red,thickness=2): AC:=line(A,C,color=red,thickness=2): AE:=line(A,E,color=red,thickness=2): BE:=line(B,E,color=red): CE:=line(C,E,color=red,thickness=2):
T:=textplot3d([[A[],"A"], [B[],"B"], [C[],"C"], [E[],"E"]], align={left,above}, font=[TIMES,18]):
S:=sphere([x0,y0,z0],R, color="Yellow"):
ABC:=polygon([A,B,C],color="LightGreen"): ABE:=polygon([A,B,E],color="LightGreen"): BCE:=polygon([B,C,E],color="LightGreen"): ACE:=polygon([A,C,E],color="LightGreen"):
display(AB,BC,AC,AE,BE,CE,T,S,ABC,ABE,BCE,ACE, scaling=constrained, axes=none, labels=["x","y","z"]);

   
Edit.              

restart;
c:=k->add(f(x[j])/mul(`if`(i<>j,x[j]-x[i],1), i=0..k), j=0..k);

# Example
c(5);

Probably this operation is simply not implemented in the  Physics[Vectors]  subpackage. But this is easy to fix using the procedure  P  as below:

restart:
with(Physics[Vectors]):
Setup(mathematicalnotation = true):

P:=proc(M,v)
add(add(M[i,j]*Component(v,j),j=1..3)*[_i,_j,_k][i],i=1..3);
end proc:

Example of use:

q_:=3*_i + 5*_j + 2*_k;
M:=<3,4,2|5,6,2|2,1,4>;
P(M,q_);

I removed unnecessary options:

restart;
F1:=1/(x + mu) + x - 1;
plots:-implicitplot(F1 = 0, mu = -8 .. 8, x = -4 .. 4, color=blue, thickness=2, gridrefine = 3);

Let a point move in the positive direction along the unit circle. It makes a complete revolution in time  T seconds and stops. Let its angular velocity decreases uniformly from the value  v0  to 0. Below we find the law of motion (the dependence of the angle of rotation  omega on the time  t) by solving the differential equation  Eq  and show the  N  positions of this point at regular time intervals and its radius-vector. Animation of this movement is also made. Of course we can change the parameters  T  and  N  as we like.

restart;
T:=10:
v:=t->-k*t+v0: # Angular velocity as a function of time t
# Law of motion: omega(t) - angle of rotation as a function of time
Eq:=diff(omega(t),t)=v(t);
dsolve({Eq,omega(0)=0});
omega:=unapply(eval(omega(t),%), t);
Sys:={v(T)=0,omega(T)=2*Pi};
solve(Sys,{v0,k});
assign(%):

with(plots): with(plottools):
N:=20:
Circle:=circle(color=blue):
t:=T*i/N:
Radius_Vector:=seq(line([0,0],[cos(omega(t)),sin(omega(t))],color=red,thickness=2), i=0..N-1):
Point:=seq(disk([cos(omega(t)),sin(omega(t))],0.025,color=red), i=0..N-1):
display(Point,Circle,Radius_Vector, axes=none);

P:=proc(t)
uses plots, plottools;
display(line([0,0],[cos(omega(t)),sin(omega(t))],color=red,thickness=2),disk([cos(omega(t)),sin(omega(t))],0.025,color=red));
end proc:

animate(P,['t'], 't'=0..10, frames=100, background=display(Circle), scaling=constrained, trace=20
);

radius_vector.mw

See help on  timelimit  command.

Denoting the elements of the matrix  T__PhPh  by  x[i,j] , i=1..3, j=1..3, we obtain a system of 3 linear equations with 9 unknowns. Obviously, in the general case, such a system has an infinite number of solutions depending on 6 parameters.

Q20210127_new.mw

restart;
M:=Matrix(3, 3, [[-1/3*a + b + c - 2*d, -1/4*a + 2/3*c, 1/15*a - 1/4*c + 2/3*d], [e + a - 2*c, 2/3*a, -1/4*a + 2/3*c], [f - 2*a, 0, 2/3*a]]);
Id:=Matrix(3,3,(i,j)->`if`(i=j,1,0));
Sys:=[seq(seq(M[i,j]=Id[i,j],j=1..3),i=1..3)];
solve(Sys);

Similarly, you can solve for any  n .

remove(t->is(op(1, t) = -1), [T]);

See help on the command  Student:-Calculus1:-ApproximateInt  for this.

Before looking for a root, it is always useful to do the simplest qualitative analysis of the equation. Since we are looking for real roots, we first find the domain of the expression  P . Obviously should be  6*lambda - 2 >= 0  or  lambda>=1/3 . So it doesn't make sense to consider  lambda < 1/3  where the expression  P  takes complex values (this is the reason for the error in your file) .

Having built the plot, we clearly see one root. This root is easily found by  fsolve  command that Carl did. Looking at the plot it seems that this root is the only one, but of course this requires a separate proof, since any plot shows the behavior of a function only in a certain finite interval. In this example, one can prove the uniqueness of the root simply by finding the derivative of  P  and checking that it is negative for all  lambda>1 . Below we also show 3 more ways to find this root.

restart;
P := (-216*sqrt(2)*lambda^2*(lambda-4/3)*arctan((1/2)*sqrt(6*lambda-2)*sqrt(2))-108*Pi*lambda^2*(lambda-4/3)*sqrt(2)+(-12*lambda+16)*(6*lambda-2)^(3/2)+(-40*lambda+32)*sqrt(6*lambda-2))/(576*lambda^2)+1/(6*lambda);
plot(P, lambda=0..2, size=[600,250]);
simplify(diff(P, lambda)); # The derivative of P
r:=fsolve(P=0); # The unique real root

# Other ways to find the root
RootFinding:-Analytic(P, re = 0.5 .. 10, im = -1 .. 1); 
Student:-Calculus1:-Roots(P, lambda = 0.5 .. 10, numeric)[]; 
Student:-NumericalAnalysis:-Bisection(P, lambda = [0.5, 10], tolerance = 10^(-9));

restart;
A:=[a,-b,c, d, -a,b, -d];
remove(t->sign(t)=-1, A);
convert(%, set);

                                 A := [a, -b, c, d, -a, b, -d]
                                            [a, c, d, b]
                                            {a, b, c, d}