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

This does  VectorCalculus[Gradient]  command.

1) In Maple  cosec  is coded as  csc .

2) The identity  6*x*cot(x) - 3*x^2*csc(x)^2=6*x*cot(x)+3*x^2*(-1-cot(x)^2)  can be proved by the direct commands  is  or  testeq :

is(6*x*cot(x) - 3*x^2*csc(x)^2=6*x*cot(x)+3*x^2*(-1-cot(x)^2));

testeq(6*x*cot(x) - 3*x^2*csc(x)^2=6*x*cot(x)+3*x^2*(-1-cot(x)^2));

                                                          true

                                                          true

restart;

f := plottools[transform]((x, y)->[x-(1/2)*y, y]):

A := plottools[polygon]([[-1, -1], [-1, 1], [1, 1], [1, -1]], style = line, color = blue):

B := POLYGONS(op(f(A)), COLOR(RGB, 1.0, 0., 0.)):

plots[display](A, B, scaling = constrained);

You have forgotten to call  plots  package. Replace  display  with  plots[display]. I have not found other errors in your code.

Exact factorization of the original polynomial is possible, but very cumbersome:

Download factorization.mw

You can use  usual polynomial interpolation.

Example:

Phi := [0, 1, 1.5, 2, 3]:

R := [1, 3, 3.5, 4.5, 6]:

f := unapply(CurveFitting[PolynomialInterpolation](Phi, R, phi), phi);

A := plot(R, Phi, coords = polar, style = point, color=blue, symbolsize=15):

B := plot(f(phi), phi = 0 .. 3, coords = polar):

plots[display](A, B, scaling=constrained);

 Addition: the same on polar grid:

plots[display](A, B, axiscoordinates = polar);

with(plots):
with(plottools):
with(SolveTools):
with(DocumentTools):
Sol := solve({(x-1)^2+(y-1/.8)^2 = (1/.8)^2, (x-1.2/(1+1.2))^2+y^2 = 1/(1+1.2)^2}, [x, y]):
x0 := rhs(Sol[2, 1]): y0 := rhs(Sol[2, 2]):
Point2 := disk([x0, y0], 0.02, color = red):
R := sqrt(x0^2+y0^2):
Do(sc = evalf(R)):
cir := circle([0, 0], sc):
display(cir, Point2, view = [-1 .. 1, -1 .. 1], scaling = constrained);

Expr:=-tau*alpha*sigma*omega*(tau-omega^(-sigma))/(s*L*(1-gamma-tau^2+tau^2*gamma))

-(sigma*alpha)/s*(1-tau*omega^sigma)/(omega^(sigma-1)*(1-gamma)*L*(1/tau-tau));

A:=factor(denom(op(1, Expr)));

subs(op(1,Expr)=numer(op(1,Expr))*1/A, Expr);

X:=[10,15,20,25,30,35,40]: 

Y:=[2.01,4.87,6.50,9.875,14.00,18.90,24.40]:

 f:=x->`if`(is(x in X), Y[ListTools[Search](x,X)], undefined):

Examples:

f(35);    f(1);

                  18.90

                undefined