Carl's code is working:

Sol:=fsolve({F,Fw,Fk,Ft}, map(`=`,{w,k,ki,T},0..infinity));
evalf(eval({F,Fw,Fk,Ft}, Sol));

@apm   Since  limit  is weaker than  MultiSeries:-limit  then before  limit  you have to use a simplifying transformation as shown by Axel Vogt. In addition, instead of  combine, you can use  expand  command:

restart;

Expr:=sqrt(-2*cos(alpha)*cos(alpha+d)+2-2*sin(alpha+d)*sin(alpha))/d:

limit(expand(Expr), d = 0, left);  limit(expand(Expr), d = 0, right);

                                                         -1

                                                          1

@Markiyan Hirnyk  Ill-conditioned and ill-posed systems of linear algebraic equations are the different things.

@Carl Love 

1) Thanks for comment on zero entries.

2)  Transpose of M does not help.  

In Maple 12:

 convert(M^%T, list)[];

@Carl Love  Done in Maple 12 Classic.

# Return the angle of a point in interval [0, 2*Pi).

Angle:= proc(P::[realcons,realcons])

local arg:= evalf(arctan(P[2],P[1]));

     `if`(arg < 0, evalf(arg+2*Pi), arg)

end proc:

SortByAngle:= (A::Matrix)->

     Matrix(sort(convert(A^%T, listlist), (P,Q)-> Angle(P) < Angle(Q)))^%T:

M:= < -1,2,3,1,1,-1 ; 0,2,3,1,0,0 >:

SortByAngle(M);

@wingjammer   My example is of a general nature, you can write in the brackets the order of columns that you need.

Specify your problem more clearly. I do not understand what you want in the general case.

@Carl Love   Your code is compact and elegant, but conceptually more complex and less effective:

t := time():  SpM(1000):  time()-t;

t := time():  Am(1000):  time()-t;

Of course, this can be done in different ways. Here are 2 other ways:

restart;

f:=x+y+z:

V:=<seq(diff(f, var), var=[x,y,z])>;

restart;

f:=x+y+z:  V:=Vector():  var:=[x,y,z]:

for i to nops(var) do

V(i):=diff(f, var[i]):

od:

V;

@Carl Love   It is sufficient to find the exact roots of the polynomial (the first command in my answer)

Another example:

factor(3*x^2-7);

factor(3*x^2-7, {sqrt(3), sqrt(7)});

 -op(1,%)*3*``(expand(op(2,%)/sqrt(3)))*``(-expand(op(3,%)/sqrt(3)));

@Carl Love   In your example we can do exact symbolic factorization:

A := factor(x^4+1, sqrt(2));

B := ``(op(1, A)*op(2, A))*op(3, A);

@Markiyan Hirnyk  My program has the following differences from these:

1) The program works for systems of all sizes, and these tutors only if the every dimension of the matrix  no greater than 5.

2) My program in the during of Gauss - Jordan elimination , simulating manual solution produces some simplifying transformations. For example, solving the system  {19*x+20*y=39, x+y=2}, to get rid from fractions, useful to swap positions of equations.

@Axel Vogt   Yes, if we increase Digits, the error disappears:

Digits := 30:

KsT := 0.10e11: KaT := 4.6*0.10e6: KsE2 := 3.14*0.10e10: KaE2 := 4.21*0.10e6: KsDHT := 3*0.10e6: KaDHT := 3.5*0.10e6:

eq1 := E2 = fE2*(1+(KaE2+Ca)/(1+KaE2*fE2+KaT*fT+KaDHT*fDHT)+KsE2*Cshbg/(1+KsE2*fE2+KsT*fT+KsDHT*fDHT)):

eq2 := T = fT*(1+KaT*Ca/(1+KaE2*fE2+KaT*fT+KaDHT*fDHT)+KsT*Cshbg/(1+KsE2*fE2+KsT*fT+KsDHT*fDHT)):

eq3 := DHT = fDHT*(1+KaDHT*Ca/(1+KaE2*fE2+KaT*fT+KaDHT*fDHT)+KsDHT*Cshbg/(1+KsE2*fE2+KsT*fT+KsDHT*fDHT)):

E2, T, DHT, Ca, Cshbg := 1, 2, 3, 4, 5:

solve({eq1, eq2, eq3}, useassumptions) assuming fE2 > 0, fT > 0, fDHT > 0;

     {fDHT = 0.300792608884171123405264656623e-5, fE2 = 3.02872370161763976538808187091*10^(-10), fT    = 1.90199449982478958012781172005*10^(-10)}

@velimir74  I do not know the cause of this extra solution. In fact, this is a bug, because easy to verify that this is not a solution:

KsT := 0.10e11: KaT := 4.6*0.10e6: KsE2 := 3.14*0.10e10: KaE2 := 4.21*0.10e6: KsDHT := 3*0.10e6: KaDHT := 3.5*0.10e6:

eq1 := E2 = fE2*(1+(KaE2+Ca)/(1+KaE2*fE2+KaT*fT+KaDHT*fDHT)+KsE2*Cshbg/(1+KsE2*fE2+KsT*fT+KsDHT*fDHT)):

eq2 := T = fT*(1+KaT*Ca/(1+KaE2*fE2+KaT*fT+KaDHT*fDHT)+KsT*Cshbg/(1+KsE2*fE2+KsT*fT+KsDHT*fDHT)):

eq3 := DHT = fDHT*(1+KaDHT*Ca/(1+KaE2*fE2+KaT*fT+KaDHT*fDHT)+KsDHT*Cshbg/(1+KsE2*fE2+KsT*fT+KsDHT*fDHT)):

E2, T, DHT, Ca, Cshbg:=1, 2, 3, 4, 5:

eval({eq1, eq2, eq3}, {fE2 = 7.15539742720000*10^10, fT = 3.30345789526199, fDHT =

4.72656319086592*10^14});

                               {1 = .7155397427e11, 2 = 3.303457895, 3 = .4726563191e15}

@wingjammer

UniqueColumns1:=proc(A::Matrix)

local L;

uses LinearAlgebra, ListTools;

L:=[Categorize((x,y)->Equal(x,y), [seq(A[..,i],i=1..ColumnDimension(A))])];

seq(L[i,1], i=1..nops(L));

end;

Your example:

Q := Matrix([[1,1,1,0,0,0,-2,-1], [0,0,0,1,1,1,0,-3]]);

UniqueColumns1(Q);

Addition: another way.

UniqueColumns2:=proc(A::Matrix)

uses LinearAlgebra, ListTools;

MakeUnique([seq(A[..,i], i=1..ColumnDimension(A))], 1, (x,y)->Equal(x,y));

end; 

@Carl Love   

sqrt(139.^2+105^2);

                      174.2010333