In real domain the two-argument function  arctan(y,x)  returns the polar angle  alpha  of any point with the coordinates  (x,y)  on the plane in the range  -Pi < alpha <= Pi . Below are all the variants as it expressed through the standard function arctan(x) :

arctan(y,x) assuming x>=0,y>=0;
arctan(y,x) assuming x>=0,y<0;
arctan(y,x) assuming x<0,y>=0;
arctan(y,x) assuming x<0,y<0;

If I understand your problem correctly, then you want to get derivatives of different orders from the original expression  u[0]  using the formulas for differentiating the sum and the product at each step of the for-loop. See the following code for your specific example:

restart;
u[0] := exp(-x^2);
n:=10:
d := r->diff(coeff(r,u[0]),x)*u[0]+coeff(r,u[0])*diff(u[0],x):

for k from 1 to n do 
if type(u[k-1],`+`) then u[k] := map(d,u[k-1]) else 
u[k] := d(u[k-1]) end if; 
print(u[k]); 
end do:

Sx:=1/sqrt(2)*Matrix([[0,1,0],[1,0,1],[0,1,0]]);
lam,v:=LinearAlgebra:-Eigenvectors(Sx);
map(r->Vector(r/~LinearAlgebra:-VectorNorm(Vector(r),2)),convert(v^%T,listlist)); 

The problem is easy to solve if you just cut out regions with asymptotes from the plot:

restart;
func:=unapply(-tan((-4*t+x+2*y)/(2*sqrt(15)))/(2*sqrt(2)), x,y,t):
y1:=solve((-4*t+x+2*y)/(2*sqrt(15))=-Pi/2,y);
y2:=solve((-4*t+x+2*y)/(2*sqrt(15))=Pi/2,y);
y3:=solve((-4*t+x+2*y)/(2*sqrt(15))=3*Pi/2,y);
y4:=solve((-4*t+x+2*y)/(2*sqrt(15))=5*Pi/2,y);
A:=plot3d(func(x,y,4), x=-10..10, y=eval(y1+0.1..y2-0.1,t=4), style=surface):
B:=plot3d(func(x,y,4), x=-10..10, y=eval(y2+0.1..y3-0.1,t=4), style=surface):
C:=plot3d(func(x,y,4), x=-10..10, y=eval(y3+0.1..y4-0.1,t=4), style=surface):
plots:-display(A,B,C, scaling=constrained);

restart;
eq := cos(2*x)-sin(x)*(sin(x)^2+1)^(1/2)+cos(x)^2*sin(x)/(sin(x)^2+1)^(1/2):
solve(subs(x=arcsin(y), eq));

By default, Maple works in the complex domain and, of 3 values for a cube root, returns the value with the smallest argument (the principal value of a root). When working in the real domain, use the  surd  command:

a:=(-8)^(1/3);
evalc(a);
surd(-8,3);
plot(surd(x,3), x=-8..8, size=[1000,250], scaling=constrained, labels=[x,surd(x,3)]);

An alternative for  surd  would be to call the package  RealDomain :

with(RealDomain):
plot(x^(1/3), x=-8..8, size=[1000,250], scaling=constrained, labels=[x, x^(1/3)]);

It is best not to use the command to draw a surface of revolution. Below, to construct the intersection, the cone and the plane are given by explicit equations in 3D, and both spheres are given by graphic primitives:

restart;
# f := x -> 0.5*x;
# small sphere
# g := x -> sqrt(0.8944271908^2 - (x - 2)^2)
# bigger sphere 
# h := x -> sqrt(3.354101965^2 - (x - 7.5)^2)

with(plots): with(plottools):
Eq:=y=1.216350358*x - 3.841106394:
Plane:=implicitplot3d(Eq, x=0..12, y=-5..5, z=-5..5, style=surface, color=green):
Cone:=plot3d([2*r,r*cos(phi),r*sin(phi)], phi=0..2*Pi, r=0..5, style=surface, transparency=0.5):
small_sphere:=sphere([2,0,0],0.8944271908, style=surface, color=gold):
bigger_sphere:=sphere([7.5,0,0],3.354101965, style=surface, color=gold):
Ellipse:=intersectplot(surface([2*r,r*cos(phi),r*sin(phi)], phi=0..2*Pi, r=0..5),surface(Eq, x=0..12, y=-5..5, z=-5..5), thickness=2): 
display(Cone,Plane,small_sphere,bigger_sphere,Ellipse, axes=normal, scaling=constrained);

Edit. You can make the cutting plane partially transparent using the option  style=line  instead of  style=surface  or the transparency  option.

Right-click on this picture, select  Export  from the pop-up menu and save the picture in the desired format.

The following code is working:

restart; 
A := Matrix([[10, 2, 0, 0, 0], [1, 4, 2, 0, 0], [0, 1, 6, 2, 0], [0, 0, 2, 7, 3], [0, 0, 0, 2, 18]], datatype = float, storage = sparse);
b := Vector([11, 20, 34, 56, 77], datatype = float);
X := LinearAlgebra:-LinearSolve(A, b, method = SparseIterative);

Edit. The help says that the matrix can be nonsymmetric.

We see that this equation has 1 real and 2 complex roots:

[solve(5*x^3-x^2+x-1=0, explicit)];
evalf~(%);

Example:

restart;
P:=randpoly([x,y]);
C:=[coeffs(P,[x,y], 't')];
T:=[t];
TC:=zip(`[]`,T,C);
L:=[seq([[degree(p[1],x),degree(p[1],y)],p[2]], p=TC)];

I think this integral can only be calculated numerically when all parameters are given. 

Example:

restart;
Cap_Unit_Length := Int((epsilon_0*epsilon_r*d*sqrt(1+(T*Pi*S)^2))/(epsilon_r*S*theta*tan(theta/2)+2*t),theta=0..Pi);
evalf(eval(Cap_Unit_Length, [epsilon_0=1,epsilon_r=2,d=3,T=4,S=5,t=6]));

You can easily do this by removing the term  O(x^6)  and adding the required term (here  xi  is between  0  and  x ):

restart;
convert(series(cos(x),x=0, 5), polynom)+(D@@6)(cos)(xi)*x^6/6!;

See the general formula in  https://en.wikipedia.org/wiki/Taylor%27s_theorem

Edit.

I think that different forms of answers are related to different methods used in the calculation. Use method=FTOC option for simplicity in this example:

int(2*Pi*x*sqrt(1+x^(-4/3)/9), x=0..1, method=FTOC);

FTOC (abbreviation) = fundamental theorem of calculus

Since the interval for the unknown  h  (in which we are looking for a solution) has already been indicated  h>=10 and h<=20, everything can be done in 1 step:

restart;
f := x -> sqrt(-x^2+20*x):
solve(2*Pi*int(f(x)*sqrt(1+diff(f(x),x)^2),x=0..h)=1005) assuming h>=10,h<=20;

                                         201/(4*Pi)