Use  factors  command for this:

f:=x*(x+1)^4*(x^2+x+1)*(x^3+x^2+1)^5:
factors(f);
map(t->op(1,t), %[2]);
sort(%);
                         

I usually use  is  command for such checks:

is(2*cos(phi)^2-1=cos(2*phi));
                                                             true

And more often I look  here , if I don’t know any rare identity.

Edit.

Maple 2018.2 returns the result in the form of finite sums and products (without hypergeometric functions):

restart;
int(sin(x)^n, x) assuming n::posint;
eval(%, n=10);
value(%);
int(sin(x)^10, x); # The direct calculation

restart;
Ke := Matrix(4, 4, {(1, 1) = (12*I)*E/l^3, (1, 2) = (6*I)*E/l^2, (1, 3) = -(12*I)*E/l^3, (1, 4) = (6*I)*E/l^2, (2, 1) = (6*I)*E/l^2, (2, 2) = (4*I)*E/l, (2, 3) = -(6*I)*E/l^2, (2, 4) = (2*I)*E/l, (3, 1) = -(12*I)*E/l^3, (3, 2) = -(6*I)*E/l^2, (3, 3) = (12*I)*E/l^3, (3, 4) = -(6*I)*E/l^2, (4, 1) = (6*I)*E/l^2, (4, 2) = (2*I)*E/l, (4, 3) = -(6*I)*E/l^2, (4, 4) = (4*I)*E/l}):
k:=E*I/l^3:
``(k).(Ke/k);

You can take any initial value (also any range) for the parameter yourself.
Example:

restart;
plots:-animate(plot, [m*x, x=-10..10, -10..10], m=1..5, frames=100);

Here is a simple procedure for this and an example:

Addition. Below is an example of another method (without a custom procedure):
 

A:=Matrix(2,3,(i,j)->LinearAlgebra:-RandomVector(2,generator=-9..9));
B:=Matrix(3,4,(i,j)->LinearAlgebra:-RandomVector(2,generator=-9..9));
subsindets(A.B, `*`,t->`.`(op(t)));

Download P.mw

At the beginning of the code, you enter a recursive function  Tm . But your method does not work. Take a look:
 

restart;
Tm:=(t,m)-> 2*t*Tm(t,m-1)-Tm(t,m-2):
Tm(t,0):=1:
Tm(t,1):=t:
Tm(3,4);  # Example

     Error, (in Tm) too many levels of recursion

You can use the following procedure to correctly determine:
 

restart;
Tm:=proc(t,m)
option remember;
if m=0 then return 1 else
if m=1 then return t else
2*t*Tm(t,m-1)-Tm(t,m-2) fi; fi;
end proc:
Tm(3,4);  # Example

                                   577

You can also get an explicit formula for TM :
 

restart;
Tm:=unapply(rsolve({Tm(t,m)=2*t*Tm(t,m-1)-Tm(t,m-2), Tm(t,0)=1, Tm(t,1)=t}, Tm(t,m)),t,m);
Tm(3,4);
simplify(%);

We can use the inert form integral for this:

G := Int( exp( -epsilon * ( x^4 + x^2 ) ), x = -10 .. 10);
plot( G, epsilon=0 .. infinity );
         

Download QuestionSimulation_new.mw

RootFinding:-Analytic(exp(x)*sin(x)-exp(1), re=0..100, im=-1..1);
plot(exp(x)*sin(x)-exp(1), x=0..100,-3..3, size=[900,300]);  # Check

Notice that the exponential function with the base  e  must be encoded in Maple as  exp(x)  (not  as  e^x ). e is just a symbol in Maple.

Unfortunately, I do not see the possibility to significantly speed up this calculation. I just rewrote your code as a procedure, slightly simplifying it. The procedure returns  true  if the matrix is  MDS  and  false  otherwise. The example of a matrix of 12 order is considered. This calculation took about 35 minutes on my computer.

Download MDS.mw

Since the system is easily reduced to a linear system, it can easily be solved symbolically, that is, absolutely exactly:

Sys:=convert({(T[1]-T[0])/(10000-T[1]+T[0]) = -2.000000000, (T[2]-T[0])/(20000-T[2]+T[0]) = 0, (T[3]-T[0])/(50000-T[3]+T[0]) = 50, .1*T[0]+.3*T[1]+.55*T[2]+0.5e-1*T[3]-5000 = 0}, fraction);
solve(Sys);

                

To solve this problem you can use  PosIntSolve  procedure from this post. This procedure allows you to find the number of non-negative (or positive) integer solutions of a linear Diophantine equation with positive coefficients. For convenience, here are the text of the procedure and 2 examples.

restart;

PosIntSolve:=proc(L::list(posint), n::nonnegint, s::symbol:=nonneg)
local N, n0;
option remember;
if s=pos then n0:=n-`+`(op(L)) else n0:=n fi;
N:=nops(L);
if N=1 then if irem(n0,L[1])=0 then return 1 else return 0 fi; fi;
add(PosIntSolve(subsop(1=NULL,L),n0-k*L[1]), k=0..floor(n0/L[1]));
end proc:

The problem: find the number of non-negative integer solutions of equations  51*x1+53*x2+50*x3=197  and   51*x1+53*x2+50*x3=2697  .

Solutions:

PosIntSolve([51, 53, 50], 197);
PosIntSolve([51, 53, 50], 2697);
                                                             0
                                                             28
 

Replace  "5"  temporarily some symbol:

applyrule(a*x+a=a*``(x+1), subs(5=a, sqrt(5*x+5+y))):
subs(a=5, %);
                             sqrt(5*(x+1)+y)

Actually no error. When you differentiate the internal function  r*sin(theta)*(r*cos(theta)-a)^(-1)  by a-variable, then the factor (-1) occurs 2 times.