The fsolve command did not take into account the constraints that assume command imposes. To solve your problem, it is necessary to choose the starting point for the iterations. After several attempts, I have found such a point:

restart;

Digits:=20:
assume(w>0): assume(T>0): assume(k>0): assume(ki>0): assume(Mn>0): assume(Ms>0): assume(Tf>0): assume(alfa>0):
Gp:=1/exp(sqrt(w*I)):
Cf:=((T*w*I+1)/(Tf*w*I+1)):
Cpi:=(k*w*I+ki)/(w*I):
L:=Cf*Cpi*Gp:
L:=evalc(L):
F:=subs({Ms=2,Tf=(k*T^(alfa)/Mn)^(1/alfa)},evalc(abs(1+L)^2-1/Ms^2)):
F:=subs({Mn=10,alfa=0.2},evalc(F)):
Fw:=diff(F,w):
Fk:=diff(F,k):
Ft:=diff(F,T):
Sol:=fsolve({F,Fw,Fk,Ft},{w=10,k=5,ki=10,T=10});
evalf(eval({F,Fw,Fk,Ft}, Sol));

I tried to solve your system with Mathematica and it showed that the system is inconsistent:

Urgoth40, 

you made ​​a few mistakes in your code: 

1 Incorrectly placed parentheses. 

2 You wrote down only the expression for the last term, but should be the sum of terms. 

3 Since the animation parameter  n  assumes only discrete values ​​with variable pitch, it is better to use  plots[display]  with option  insequence=true  instead of  plots[animate] .

Here is a solution of your problem (12 frames allocated to every plot - approximately 1 sec):

L := [seq(5*2^i, i = 0 .. 7)];

S:=seq(plot(add(2/((2*k-1)*Pi)*sin((2*k-1)*Pi*x),k=1..n), x=-1..1) $ 12, n=L):

plots[display](S, insequence=true, scaling=constrained);

                                          L := [5, 10, 20, 40, 80, 160, 320, 640]

m:=1:  N:=2:

Matrix(N*(m+1), (i,j)->`if`(j=m+1+i, 1, 0));

k:=3:  m:=1:  P:=Matrix(m+1, [1,2,3,4]):  P1:=Matrix(m+1, fill=1)/k/(m+1):  P0:=Matrix(m+1):

Matrix([seq([P0 $ i-1, P/k, P1 $ k-i], i=1..k)]);

All the formulas you can see in the link  http://orion.math.iastate.edu/alex/166H/polar_lines_tangents.pdf

r:=theta->3-3*cos(theta):

theta1:=3*Pi/4:

a:=polarplot(r(theta),theta=0..2*Pi, color=blue, thickness=2):

phi:=-arcsin(D(r)(theta1)/sqrt(r(theta1)^2+D(r)(theta1)^2)):

d:=r(theta1)^2/sqrt(r(theta1)^2+D(r)(theta1)^2):

theta0:=theta1+phi:

b:=polarplot(d/cos(theta-theta0), theta=Pi/3..23*Pi/24, color=red):

c:=plottools[disk]([r(theta1)*cos(theta1),r(theta1)*sin(theta1)],0.1, color=red):

plots[display](a,b,c, scaling=constrained);

Carl, your code does not work for earlier versions (for example in Maple 12). Here are 2 variants for earlier versions:

M:=Matrix(3,[[1,1,1],[2,2,2],[3,3,3]]);

{parse(cat(op(map(rhs, convert(op(2,M), list)))))};

{parse(cat(op(ListTools[Flatten](convert(M, listlist)))))};

restart;

A:=Matrix( [ [1,2,3],[4,5,-1]] );

for i to 2 do

for j to 3 do

if A[i,j]>=3 then L[i,j]:=[i,j]  fi:

od: od:

convert(L, list);

 Another way:

restart;

A:=Matrix( [ [1,2,3],[4,5,-1]] );

select(x->A[x[1],x[2]]>=3, [seq(seq([i,j], j=1..3), i=1..2)]);

Procrdure  f  produces all vectors:

f:=proc(m::nonnegint,Q::Matrix)

local k, i;

global e;

for k from 0 to m do

for i from 0 to m do

e[k,i]:=binomial(m,i)/(2*m+k+1)*Q^(-1).Vector([seq(binomial(m,j)/binomial(2*m+k,i+k+j), j=0..m)]);

od; od;

end proc:

Example:

f(2, <1,3,5; 2,5,7; 3,-5,4>):

e[0,2], e[1,1];

Now I understood the meaning of the problem. Carl's solution is good, but it does not give an unequivocal result, if two different points have the same angle. Here's an easy modification of Carl's solution, which, if equal angles, sorts in ascending order of the absolute value:

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

`if`(Angle(P)=Angle(Q),evalf(P[1]^2+P[2]^2)<evalf(Q[1]^2+Q[2]^2), Angle(P) < Angle(Q))))^%T:

Example:

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

 SortByAngle1(M);

R:=<Q[.., [2, 1, 4, 3]]>;

SpM:=proc(m::nonnegint)

Matrix([seq([0 $ i, seq((-1)^k*binomial(m,i)*binomial(m-i,k), k=0..m-i)], i=0..m)]);

end:

Example:

SpM(7);

y := (r, t) -> 45*r*t;

seq(y(5, t), t=2..150);

                                                                   y := (r, t) -> 45*r*t

450, 675, 900, 1125, 1350, 1575, 1800, 2025, 2250, 2475, 2700, 2925, 3150, 3375, 3600, 3825, 4050, 4275, 4500, 4725, 4950, 5175, 5400, 5625, 5850, 6075, 6300, 6525, 6750, 6975, 7200, 7425, 7650, 7875, 8100, 8325, 8550, 8775, 9000, 9225, 9450, 9675, 9900, 10125, 10350, 10575, 10800, 11025, 11250, 11475, 11700, 11925, 12150, 12375, 12600, 12825, 13050, 13275, 13500, 13725, 13950, 14175, 14400, 14625, 14850, 15075, 15300, 15525, 15750, 15975, 16200, 16425, 16650, 16875, 17100, 17325, 17550, 17775, 18000, 18225, 18450, 18675, 18900, 19125, 19350, 19575, 19800, 20025, 20250, 20475, 20700, 20925, 21150, 21375, 21600, 21825, 22050, 22275, 22500, 22725, 22950, 23175, 23400, 23625, 23850, 24075, 24300, 24525, 24750, 24975, 25200, 25425, 25650, 25875, 26100, 26325, 26550, 26775, 27000, 27225, 27450, 27675, 27900, 28125, 28350, 28575, 28800, 29025, 29250, 29475, 29700, 29925, 30150, 30375, 30600, 30825, 31050, 31275, 31500, 31725, 31950, 32175, 32400, 32625, 32850, 33075, 33300, 33525, 33750

Your equations can only be solved numerically for the specified parameters.

Here is the example of the solution of the first equation:

restart;

b, m, y1, y2:=1, 2, 3, 4:  # Specification of the parameters

Sol:=solve({(b+m*yc)^3*yc^3/((b+2*m*yc)*(b+m*y1)*y1)-(b+m*yc)^3*yc^3/((b+2*m*yc)*(b+m*y2)*y2) = y2^2*((1/2)*b+(1/3)*m*y2)-y1^2*((1/2)*b+(1/3)*m*y1)}, {yc}):

op(map(rhs@op@fsolve, [Sol], yc, complex));  # All solutions

3.477366262,  0.9018399724+3.525754065*I,  -3.265522931+2.179059836*I,  -0.2500003440,  -3.265522931-2.179059836*I,  0.9018399724-3.525754065*I

Example:

L:=[$ 1..25];

A:=Matrix(5, L);