Will this be what you want:

restart;
S__x := sin(r__a):
S__y := cos(r__a):
S__z := sin(r__a):
lower_limit := 0: upper_limit := 10: step_size := .3: num_steps := ceil((upper_limit-lower_limit)/step_size): r_values := [seq(i*step_size+lower_limit, i = 0 .. num_steps)]:
                      
arrows := [seq(plots:-arrow([0, 0, r__a], [S__x(r__a), S__y(r__a), S__z(r__a)], color = blue), r__a in r_values)]:

plots:-display(arrows, axes = normal);

restart;
phi := proc(k,x,L)
  if (type(k,even)) then sqrt(2)*sin(Pi*k*x/L)/sqrt(L)
  else sqrt(2)*cos(Pi*k*x/L)/sqrt(L)
  end if;
end proc:
 	
A:=Int(phi(m,x,L)*_h^2/m2*diff(phi(n,x,L),x,x),x=-L/2..L/2);
value(A) assuming m::posint,n::posint,m=n;
value(A) assuming m::posint,n::posint,m<>n;

We place point  A  at the origin of coordinates, point  B  on the axis  Ox  to the right of  A , points C and  D  lie in the upper half-plane. We introduce 2 unknowns: s - the angle between the axis  Ox  and AD, t - the angle between the axis  Ox  and  BC. We use vectors for calculations. We find s  and  t  numerically by solving a system of 2 equations with 2 unknowns:

restart;
with(plottools): with(plots): 
local D:
Dist:=(X,Y)->sqrt((X[1]-Y[1])^2+(X[2]-Y[2])^2):
Angle:=(U,V)->arccos((U.V)/sqrt(U.U)/sqrt(V.V)):
A:=<0,0>: B:=<170,0>: C:=<170+170*cos(t),170*sin(t)>: D:=<180*cos(s),180*sin(s)>:
Sys:={Dist(C,D)=290,Dist(A,C)=Dist(B,D)};
Sol:=fsolve(Sys, {s=0..Pi,t=0..Pi});
C:=eval(C, Sol): D:=eval(D, Sol):
P:=curve(convert~([A,B,C,D,A],list), color=red, thickness=2):
T:=textplot([[convert(A,list)[],"A"],[convert(B,list)[],"B"],[convert(C,list)[],"C"],[convert(D,list)[],"D"]], font=[times,bold,18], align={left, above}):
display(P,T, scaling=constrained);
DA:=A-D: DC:=C-D: CD:=-DC: CB:=B-C:
eval(s,Sol); # Angle A
eval(Pi-t,Sol); # Angle B
Angle(DA,DC);  # Angle D
Angle(CD,CB);  # Angle C

If a mapping of finite-dimensional spaces is given, then the total derivative is understood as the Jacobian, represented by a matrix consisting of partial derivatives. In the case of a real-valued function, this will be a row matrix. The total differential is obtained by multiplying this matrix by the column vector of differentials of the independent variables.

restart;
# The case of a real-valued function of two variables
M:=Student:-MultivariateCalculus:-Jacobian([f(x,y)],[x,y]);
(M.<dx,dy>)[1];

u:=t(x,y)-s(x,y): # Your example
M:=Student:-MultivariateCalculus:-Jacobian([u],[x,y]);
(M.<dx,dy>)[1];

I don't know any other way than using the combination of  algsubs  with expand  (if needed).

Example:

restart;
Expr:=cos(x)^4:
expand(algsubs(cos(x)^2=1-sin(x)^2, Expr));

Use the  evalc  command for this. Unfortunately, the result is very cumbersome, because your equation contains many parameters:

restart;

expand((a+I*b)*(a-I*b));
(a+I*b)+(a-I*b);
# z = x *Iy    I need to seperate out solution to [x,y] & [x,-y]
#https://math.stackexchange.com/questions/44391/foci-of-a-general-conic-equation
eq := T*z^2 - (R + S*I)*z +G-K + H*I;
sol:=solve(eq,z);
sol1,sol2:=evalc~([%])[];
expand(sol1*sol2)=a^2+b^2;
simplify(sol1+sol2=2*a);
simplify(sol1-sol2=2*b);

Edited.

There is no need for a triple integral here. Let us consider the section of this torus by the plane  y=0 . We see that the torus is obtained by rotating the circle  (x - 4)^2 + z^2 = 4  around the axis  Oz . Using the well-known formula for the volume of a body of rotation, we get the answer  V=32*Pi^2 :

restart;
Eq:=(x^2+y^2+z^2+12)^2=64*(x^2+y^2):
plots:-implicitplot3d(Eq, x=-6.7..6.7, y=-6.7..6.7,z=-2.7..2.7, style=surface, axes=normal, numpoints=100000, scaling=constrained);
factor((lhs-rhs)(eval(Eq,y=0)));
Student:-Precalculus:-CompleteSquare(op(1,%),[x,y]);
f1,f2:=solve(%, x);
V=simplify(Pi*int(f1^2,z=-2..2)-Pi*int(f2^2,z=-2..2));

restart;
expr:= 4*cos(Pi*m)/(Pi*(2*m + 1)); 
simplify(%) assuming m::integer;
``;
4*cos(Pi*m+a)/(Pi*(2*m + 1));
expand(%) assuming m::integer;
``;
4*cos(Pi*m+2*alpha*Pi)/(Pi*(2*m + 1));
expand(%) assuming integer;

Explanation. If alpha is an arbitrary real number, then  2*alpha*Pi  will also be an arbitrary real number. Therefore, I immediately replaced this in the code with the symbol  a , then to simplify we will not need additional commands. If  alpha  is an integer, then due to the periodicity of the cosine, the result coincides with example 1. 

Code in Maple 2018.

Edited.

restart;
h:=(t,x)->(m*t^2 + 6*t - 2*x)^2/(36*g*t^2):
t:=<1,5,10,15,20>:
V:=h~(t,x);
V[1], V[3];

This equation can only be solved numerically using  the  fsolve  command. To get roots other than 0, you can first plot a graph and then indicate the boundaries of the interval in which the root you need lies:

restart;
Eq:=(5*x*cos(x^2/2)-5*sin(5*x/2)*cos(5*x/2)=0);
plot(lhs(Eq), x=0..5);
fsolve(Eq, x=1..2);

To solve this problem it is convenient to use the differentiation operator  D , as well as the well-known tangent equation  y=f(x0) + f'(x0)*(x - x0)  at the point  x=x0 :

restart;
f:=x->(1-k*x)/(1+x^2):
T:=f(3)+D(f)(3)*(x-3); # Equation of the tangent line at x=3
T:=collect(T, x);  

# Visualization for k=1
k:=1:
plots:-display(plot([f(x),T], x=-0..10, -0.3..0.1, color=[blue,red]), plots:-pointplot([3,f(3)], color=red, symbol=solidcircle, symbolsize=12));

restart;
identify(fsolve(x!=3628800, x=0..infinity));
ifactor(65536);

It is easy to check the truth of the identity   

So we get

restart;
is((2*x^2022+1)/(x^2023+x)=1/x+x^2021/(x^2022+1));
int(1/x+x^2021/(x^2022+1), x);

You can use the technique of stitching sequential animations from the post
https://mapleprimes.com/posts/207840-Combinations-Of-Multiple-Animations

restart;
Frame:=proc(t,Range,Color)
local SinArc, p, l1, l2, l3, l4, l, Circle;
uses plots, plottools;
SinArc:=plot(sin(x),x=Range[1]..t, color=Color, thickness=2);
p:=pointplot([[-2,0],[-2,sin(t)],[cos(t)-2,sin(t)],[t,sin(t)]], symbolsize=11, symbol=solidcircle);
l:=line([-2,-1.4],[-2,1.4], thickness=0);
l1:=line([-2,0],[-2,sin(t)], color=Color, thickness=2);
l2:=line([-2,0],[cos(t)-2,sin(t)]);
l3:=line([cos(t)-2,sin(t)],[t,sin(t)]);
l4:=line([-2,sin(t)],[cos(t)-2,sin(t)]);
Circle:=circle([-2,0], color=Color, thickness=2);
display(SinArc,p,l,l1,l2,l3,l4,Circle, scaling=constrained);
end proc:

Example of use:

with(plots):
A:=animate(Frame,[t,[0,2*Pi],red], t=0..2*Pi, frames=90, size=[800,200], view=-1.4..1.4):
B:=animate(Frame,[t,[2*Pi,4*Pi],blue], t=2*Pi..4*Pi, frames=90, size=[800,200], view=-1.4..1.4):
C:=animate(Frame,[t,[4*Pi,6*Pi],green], t=4*Pi..6*Pi, frames=90, size=[800,200], view=-1.4..1.4):
display([A, display(op([1,-1,1],A),B), display(op([1,-1,1],A),op([1,-1,1],B),C)], insequence, view=[-3..19,-1.4..1.4]);

restart;
local D, O;
with(Student:-MultivariateCalculus):
A := [0, 0, 0]:
B := [a, 0, 0]:
C := [a, b, 0]:
D := [0, b, 0]:
S := [0, 0, h]:
O := [x, y, z]:
lineSC := Line(S, C);
lineSD := Line(S, D);
H := Projection(A, lineSC);
K := Projection(A, lineSD);
OH:=H-O; OK:=K-O; OC:=C-O;
M:=Matrix([OH, OK, OC]);
O:=eval(O,%): # The center of the desired circle
R:=simplify(Distance(O, H));  # The radius of the circle

Addition - visualization made from calculation results:

with(plots): with(plottools):
a:=2: b:=3: h:=4:
BAD:=curve([B,A,D],linestyle=3,color=blue):
BCD:=curve([B,C,D],thickness=2,color=blue):
SA:=line(S,A,linestyle=3,color=blue, thickness=2): SB:=line(S,B,color=blue, thickness=2): SC:=line(S,C,color=blue, thickness=2): SD:=line(S,D,color=blue, thickness=2):
AH:=line(A,H,linestyle=3,color=black): AK:=line(A,K,linestyle=3,color=black):
p:=Plane(C,S,D):
s:=(x-O[1])^2+(y-O[2])^2+(z-O[3])^2=R^2:
c:=intersectplot(GetRepresentation(p),s, x=-1..a+0.7, y=-1..b+0.7, z=-1..h+0.7, color=red, thickness=3):
T:=textplot3d([[A[],"A"],[B[],"B"],[C[],"C"],[D[],"D"],[S[],"S"],[H[],"H"],[K[],"K"]], font=[times,roman,18], align=left):
display(BAD,BCD,SA,SB,SC,SD,AH,AK,c,T, axes=none, scaling=constrained);