You can get several more solutions from one solution if you use transformations that transform a cubic integer lattice into itself. They will be, for example, reflections about the coordinate planes, rotations around the coordinate axes at angles that are multiples of 90 degrees, and some others. But in terms of geometry, you end up with the same pair of vectors. New solutions can be easily obtained using the brute force method. First we find the list L and then from the list  L, the list  L1  is obtained, where none of the coordinates is equal to 0:

restart;
k:=0:
for a from 0 to 9 do
for b from a to 9 do
for c from b to 9 do
for x from -9 to 9 do
for y from -9 to 9 do
for z from -9 to 9 do
if 4*(a*x+b*y+c*z)^2=3*(a^2+b^2+c^2)*(x^2+y^2+z^2) and a^2+b^2+c^2<>0 and x^2+y^2+z^2<>0 and a*x+b*y+c*z>0 then k:=k+1;
L[k]:=[<a,b,c>,<x,y,z>] fi;
od: od: od: od: od: od:
L:=convert(L,list):  
nops(L);
L1:=select(t->`and`(seq(t[1][i]<>0,i=1..3),seq(t[2][i]<>0,i=1..3)),L):
nops(L1);

The very first pair of vectors in the lists  L  and  L1  is your original example, in which the first and third coordinates are reversed (reflection of space relative to the plane  x=z ). To reduce the search, I chose the vector  <a,b,c>  with non-decreasing coordinates.

This can be done in many ways. Here's one:

A := [1,0,3,4,5,5]: B := [0,1,2,6,3,6]:
nops(select(`>`, A - B, 0));

                                           3

The simple recursive procedure does the work:

restart;
a:=proc(n)
local m;
option remember;
if n=1 then return 1 elif n=2 then return 0 else
for m from 1 to n-2 do
if a(n-1)=a(n-1-m) then return m fi; od;
min(seq(abs(a(n-1)-a(k)),k=1..n-2)) fi;
end proc:

Examples of use:

seq(a(n), n=1..100);
plots:-listplot([%], style=point, symbol=solidcircle, color=red);

Edit. The procedure above works about 2 times faster than tomleslie's method.

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