eval(x[1,2,1], Sol[2]);

                                                               1

I have assumed that the angle is 35 degrees (by default  in Maple angles in radians). The system has 2 solutions in the range  -Pi<=x<=Pi :

a:=35*Pi/180:
Eq1:=-sin(a)*T1+T2=0:
Eq2:=40-T1*cos(a)=0:
Eq3:=T3*sin(x)-T2=0:
Eq4:=50+cos(x)*T3=0:
Sol:=solve({Eq1,Eq2,Eq3,Eq4}, explicit);  # Symbolic (exact) answer
evalf(Sol);  # Numerical (approximate) answer

Edit. We can slightly simplify the resulting symbolic answer  Sol  by :

simplify([Sol]);
 

I made a few corrections in the latter group.

quest1.mw

restart;
f:=y^3+x^2:
M:=maximize(f, x=0..1, y=-1..1);
A:=plot3d(f-M, x = 0 .. 1, y = -1 .. 1, style=surface, color=khaki, filled):
F:=plottools:-transform((x,y,z)->[x, y, z+M]):
plots:-display(F(A),  orientation=[45,75]);

List:=[13,16,16,29,34,33,33,12,22,26,25,25,25,11]:
ListTools[Search](25, List);  # The position of the first number 25  in the list  List
List:=subsop(%=NULL, List);

                                              11
   [13, 16, 16, 29, 34, 33, 33, 12, 22, 26, 25, 25, 11]
 

Here is another way similar Preben's one, but shorter and without tricky commands as freeze, evalindets and etc.

restart;
sys:={u+v = a , u-v = b};
solve(sys, {u,v});
eval(%, {a=<-2,4>, b=<-3,6>});

The problem is easily reduced to solving a system of four equations with four unknowns:

u:=<u1,u2>:  v:=<v1,v2>:
Eq1:=u+v =~ <-2,4>;  Eq2:=u-v = ~<-3,6>;
solve({Eq1[1], Eq1[2], Eq2[1], Eq2[2]});
eval([u, v], %)[ ];

Addition. Here is a version with  Equate  command:

u:=<u1,u2>: v:=<v1,v2>:
Sys:=op~({Equate(u+v,<-2,4>), Equate(u-v,<-3,6>)});
solve(Sys);
u, v:=eval([u,v], %)[ ];

This is your error rather than Maple. You have entered a non-linear equation, and Maple only solves (step by step) linear equations. See Instructions below in the window of this application.

Addition.  See this application  http://www.universalmathsolver.com/    especially here

http://download.cnet.com/Free-Universal-Algebra-Equation-Solver/3000-2053_4-75211972.html?part=dl-&subj=dl&tag=button

I did not understand that your procedure solveEeaMatrix should be doing.

Here is an implementation of the extended Euclidean algorithm with Maple. The procedure  ExtendedEuclid  returns the greatest common divisor  d  of numbers  a  and  b  and  2 numbers  x  and  y  that  a*x+b*y=d

ExtendedEuclid:=proc(a::nonnegint, b::nonnegint)
local d, x, y, x1, y1, x2, y2, a1, b1, q, r;
if a<b then error "Should be a>=b" fi;
if b=0 then d:=a; x:=1; y:=0; return [d,x,y] fi;
x2:=1; x1:=0; y2:=0; y1:=1; a1:=a; b1:=b;
while b1>0 do
q:=floor(a1/b1); r:=a1-q*b1; x:=x2-q*x1; y:=y2-q*y1;
a1:=b1; b1:=r; x2:=x1; x1:=x; y2:=y1; y1:=y;
od;
d:=a1; x:=x2; y:=y2;
[d,x,y];
end proc:
 

Example of use:

ExtendedEuclid(30,12);

                                                          [6, 1, -2]

We see that  30*1 + 12*(-2) = 6

If four points lie in the same plane, then such the unique sphere does not exist. This case should be provided for the procedure:

getEq := proc(L::listlist) 
local p, S;
uses geom3d; 
seq(point(p || j, L[j]), j = 1 .. 4);
if AreCoplanar(p1,p2,p3,p4) then return `Not exist` else 
Equation(sphere(S, [p1, p2, p3, p4], [x, y, z])) fi; 
end proc:

Example of use:

map(getEq, L);

Edit.

If  A  and  B  are 2 points then  t*A+(1-t)*B  ( 0<=t<=1 )  is the segment  AB .

Example:

A:=[1,1]:  B:=[3,2]:
plot([op(t*~A+(1-t)*~B), t=0..1], color=red, thickness=3, view=[0..4, 0..3], scaling=constrained);

Addition.  You can also use  geometry:-segment  command or  plottools:-line  command.

seq(point(M||i, L[i]), i=1..nops(L));
seq(Equation(plane(Q||i, [M||i, n], [x,y,z])), i=1..nops(L));

1. Should be  coordinates  instead of   cordinates .

2. It will be easier to find the intersection of two curves by  solve  command:

solve([x^2+y^2-6600*x-4400*y+15730000=12100, x^2+y^2-6820*x-4840*y+17484500=6400], explicit);

We see that these circles do not intersect.

Visualization:

plots:-implicitplot([x^2+y^2-6600*x-4400*y+15730000=12100, x^2+y^2-6820*x-4840*y+17484500=6400], x=3100..3600,y=2000..2500, scaling=constrained, gridrefine=3);

The procedure  FP  solves your problem for a plane  P . Using  map  command you can solve the problem for the list of several planes.

FP:=proc(P::`=`, var::list, A::list, B::list)
local n1, n2, n, k;
uses LinearAlgebra;
n1:=convert(map2(coeff, lhs(P), var), Vector);
n2:=convert(B-A, Vector);
n:=n1 &x n2;
k:=sort(primpart(n.convert(var-A, Vector)));
k*signum(lcoeff(k)) = 0;
end proc:

Examples of use ( L  is your list above):

FP(L[5], [x,y,z], [1,5,3], [-2,-4,3]);

map(FP, L, [x,y,z], [1,5,3], [-2,-4,3]);

Addition. The procedure  FP  will not give the correct answer if the plane  P is perpendicular to the vector AB, because in this case, there is no an unique solution. You are free to supplement the procedure for this case.

If we have  a*x^2+b*x+c=0 and  a>0 and  b>0 and  c>0  then if  d=b^2-4*a*c  there are 3 cases only:

1. If  d>0  then  both roots are different and negative.

1. If  d=0  then  both roots (actually 1 root)  are equal and negative.

1. If  d<0  then  there are no real roots.

All this is easy to prove, for example using the Vieta formula.

Visualization of all the cases:

A:=plot([2*x^2+4*x+1, x^2+4*x+4, x^2+3*x+6], x=-4..1, color=[red,blue,green, thickness=2]):
T:=plots:-textplot([[-4,4.5,"d = 0"],[-4,10.5,"d < 0"],[-3.5,16,"d > 0"]], font=[roman,bold,16]):
plots:-display(A,T);