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

It is well known that a square system has a unique solution if and only if it's main determinant is not 0.

restart;

A:= < k+3, 2, k-4, 3;

      0,   2,  -9, 5;

      0,   0, k^2+k-2, k-1 >;

B:=A[..,1..3];

V:=A[..,4];

LinearAlgebra[Determinant](B);

Sol:=[solve(%)];

We have a unique solution for any number  k  not included in this list.

Next solve separately for each number from  list  Sol:

for i from 1 to 3 do

'k'=Sol[i],  LinearAlgebra[LinearSolve](op(eval([B,V],'k'=Sol[i])), free=C);

od;

Thus for  k=-3  and  k=1  we have infinitely many solutions, for  k=-2  no solutions

You made ​​a few syntax errors, which I corrected.
Your nonlinear system has 5 parameters. Maple can not solve it symbolically for all values ​​of the parameters, but Maple can  solve this system for the specified parameters.

restart;

f := proc (E2, T, DHT, Ca, Cshbg)

local KsT, KaT, KsE2, KaE2, KsDHT, KaDHT, eq1, eq2, eq3;

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

solve({eq1, eq2, eq3})

end proc:

Example:

f(1, 2, 3, 4, 5);

We got 10 solutions.

Procedure  UniqueColumns  returns  the unique columns of a Matrix.

UniqueColumns:=proc(A::Matrix)

local L;

uses LinearAlgebra, ListTools;

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

op(map(op,select(x->is(nops(x)=1), L)))

end;

Your example:

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

UniqueColumns(Q);

Write your function as a procedure:

f:=proc()

your expression

end;

Example:

f:=proc()

add(args[i]^2, i=1..nargs)

end;

f(1, 2, 3);

       14