Use  simplify  with side relation for this:

restart;
simplify(exp(-a), {exp(a)=p});

Also

f:=exp(-a)*(b+c)*1/k;
simplify(f, {exp(a)*k=p});

The calculations on my machine (Maple 2018) took only 23 seconds. Found 64 solutions:

restart;
ts:=time():
k:=0:
for x from 3 to 30 do
for y from 2 to x-1 do
for z from 1 to 30 do
for a from 3 to 10 do
for b from 2 to a-1 do
for n from 2 to 10 do
for m from 2 to n do
if igcd(m,n)=1 and igcd(a,b,m)=1 then 
if 1/x+1/y=1/z and a/x+b/y=m/n then k:=k+1; L[k]:=[x,y,z,a,b,m,n]
fi; fi;
od: od: od: od: od: od: od:
L:=convert(L, list);
nops(L);
time()-ts;

Equality  x+y=1  is possible only if  y=1-x . But then the contradiction disappears:

diff(x+(1-x)=1, x);

                          0 = 0

For other values of  y , there is no any equality. Therefore, it is not surprising that the derivatives do not coincide.

General conclusion: if some equality  f(x)=g(x)  is an identity (that is, it is true for all values of the variable  x) then  f(x)=g(x) => diff(f(x),x)= diff(g(x),x) . If equality is not an identity, then both cases are possible  diff(f(x),x)= diff(g(x),x)  as well as diff(f(x),x)<>diff(g(x),x) .

series(1/(1-x^2), x=1);

Or

numapprox:-laurent(1/(1-x^2), x = 1)

This expansion is hold in the ring  0<|x-1|<2  containing no singular points  x=1  and  x=-1 .

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 .