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

If we assume that a solution exists and substitute into the equation  x = 0, we obtain  theta (0) = 0. But then, solving with these initial conditions, we have only the trivial solution:

eq:=(1+theta(x))*diff(theta(x),x)+x*diff(theta(x),x)^2+x*(1+theta(x))*diff(theta(x),x,x)-theta(x)^3-diff(theta(x),x)=0;

ics:=theta(0)=0, D(theta)(0) = 0;

Sol:=dsolve({eq, ics}, numeric);

plots[odeplot](Sol,[x,theta(x)], x=0..1, axes=box);

Here is step by step solution. I used my procedure  JordanGausse . This procedure step by step solves any system of linear equations by Gauss - Jordan elimination. In the procedure, the Russian text is used. After I prepare the English version, I will submit  code of the procedure to the mapleprimes forum.

Unfortunately, when loading the Russian text was distorted. Here is a translation into English:

Step 1: write the system in matrix form

Step 2: divide row 1 by 4

Step 3: row 1, multiplied by -2, add to row 2

Step 4: row 1, multiplied by -4, add to row 3

Step 5: swap the positions of columns 2 and 3

Step 6: divide row 2 by 3/2

and so on ...

Download Gauss-Jordan.mws