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

Here is a more efficient way:

restart;
GetParts:=proc(L::list)
local L1, n, a, S, L2, L3, L4;
 uses combinat, ListTools;
 L1:=sort(L); n:=nops(L); a:=`+`(L[])/2;
 S:={seq([i,L1[i]], i=1..n)};
 L2:=[seq(op(choose(S,j)), j=1..floor(n/2))];
 L3:=select(s->`+`(seq(s[i,2], i=1..nops(s)))=a, L2);
 L4:=map(s->[s,S minus s], L3);
 {map(t->[[seq(t[1,i,2], i=1..nops(t[1]))],[seq(t[2,i,2], i=1..nops(t[2]))]], L4)[]};
 map(t->t[1],[Categorize((x,y)->convert(x,set)=convert(y,set),convert(%, list))]);
 end proc:

Examples of use:

CodeTools:-Usage(GetParts([3,1,1,2,2,1]));
CodeTools:-Usage(GetParts([3,1,1,2,2,1,4,4,8,8]));
CodeTools:-Usage(GetParts([3,1,1,2,2,1,4,4,8,8,10,10,12,12,12]));
 

Edit.

The procedure getCoeff returnes the coefficient in front of the monomial  t  in the multivariate polynomial  P :

getCoeff:=proc(P, T, t)
local L, H, i, k:
L:=[coeffs(P, T, 'h')]: H:=[h]: k:=0:
for i from 1 to nops(H) do
if H[i]=t then k:=L[i] fi:
od:
k;
end proc:

Examples of use:

f := a^2*b^2*c^2 + 2*a^2*b^2 + 2*a^2*c^2 + 2*b^2*c^2 + a^2 - 6*a*b - 6*a*c + b^2 - 6*b*c + c^2 + 8:
getCoeff(f, [a,b,c], a^2*b^2*c^2);
getCoeff(f, [a,b,c], a^2*b^2);
getCoeff(f, [a,b,c], a*b);
getCoeff(f, [a,b,c], a^2);
getCoeff(f, [a,b,c], 1);

                                               1

                                               2

                                              -6

                                               1

                                               8

Edit.

1. If a matrix has one column, then it is better to set it as a vector not as a matrix.

2. For use at different  n  it is more comfortable to define  b as a procedure:

b := n -> Vector(n,  i->add(1/(i+j-1), j=1..n)):
 

Example of use:

seq(b(n), n=1..7);

It is better to avoid the construction  L:=[op(L), ...] , because it is very inefficient. The next version is faster and works properly:

restart; 
n := 0:
for a from -10 to 10 do
for b from -10 to 10 do
for c from -10 to 10 do
if a*b*c <> 0 then n:=n+1; k := (-a*b*c+a*b*z+a*c*y+b*c*x)/igcd(a*b, b*c, c*a, a*b*c); L[n] := [[[a, 0, 0], [0, b, 0], [0, 0, c]], k*signum(lcoeff(k)) = 0] fi;
od: od: od: 
L:=convert(L, list);
nops(L);