If you don't like an expression with sine instead of cosine, you can easily convert sine to cosine using the  convert  command:

restart;
eqaux_2 := l__bb = L__aa0 + L__aa2 * cos(2*theta - 4/3*Pi);
eqaux_2 :=convert(eqaux_2, cos);

Here we get the expression with a cosine even simpler than the original one  ( cos(2*theta - 4/3*Pi) = cos(2*theta+2*Pi/3)  is an identity).

You can use the  plots:-display  command to place multiple graphs on one plot.

Example:

restart;
sys := {diff(x(t), t) = y(t), diff(y(t), t) = -x(t)-(1/2)*y(t)}:
A:=DEtools:-DEplot(sys, [x(t), y(t)], t = 0 .. 15, [[x(0) = 1, y(0) = 0]]):
B:=plot(-x^2, x=-0.5..1):
plots:-display(A,B);

Seems like a bug. But  assume  on a separate line works even without  simplify  and  real  (in Maple 2018.2):

restart: with(Physics):
assume(a>b[2]):
csgn(a-b[2]);  # OK

The code below finds some examples of tetrahedra with these properties. Your example is the first on the list  L :

restart:
CenterOfInSphere:=proc(P1::Vector,P2::Vector,P3::Vector,P4::Vector)
local s1, s2, s3, s4, s;
uses LinearAlgebra;
s1:=1/2*Norm(CrossProduct(P3-P2,P4-P2),2, conjugate=false):
s2:=1/2*Norm(CrossProduct(P3-P1,P4-P1),2, conjugate=false):
s3:=1/2*Norm(CrossProduct(P2-P1,P4-P1),2, conjugate=false):
s4:=1/2*Norm(CrossProduct(P2-P1,P3-P1),2, conjugate=false):
s := s1+s2+s3+s4:
simplify(1/s*(s1*P1+s2*P2+s3*P3+s4*P4));
convert(%, list);
end proc:

N:=20:
k:=0:
for a from 1 to N do
for b from a to N do
for c from b to N do
C:=CenterOfInSphere(<0, 0, 0>, <a, 0, 0>, <0, b, 0>, <0, 0, c>);
R:=C[1];
if R::integer then k:=k+1; L[k]:=[[a,b,c],R] fi;
od: od: od:
L:=convert(L, list);

L := [[[3, 5, 14], 1], [[3, 6, 9], 1], [[4, 4, 8], 1], [[6, 12, 18], 2], [[6, 13, 16], 2], [[8, 8, 16], 2], [[8, 10, 11], 2], [[12, 14, 18], 3]]

You can increase the  N  number to get more solutions.

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);