If you have specific points, you can use  geom3d  package.

Example:

restart;
with(geom3d):
point(P, 5, 6, 7);
point(A, 1, -2, 3);
point(B, 2, 4, 1);
line(AB, [A,B]);
projection(Q, P, AB);
detail(Q);
                               

A:=isolve(n+k=2*c);
is(A<>NULL and eval(n, A)::name);
                          A := {c = _Z1, k = 2*_Z1-_Z2, n = _Z2}
                                                       true

Above  _Z1  and  _Z2  are any integers.

You can use  PolyhedralSets  package to get the vertices and edges of your polygon. This also works in the multidimensional case.

Your example (failed to load the contents of the file, so the image is shown, the link to the file below):

You can extract all vertices in one list with  Relations  command:
V:=Vertices(ps):
Relations~(V);
map(v->eval([x,y],v), %);
                  [[y = -1, x = -1], [y = 1, x = -1], [y = -1, x = 1], [y = 1, x = 1]]
                                       [[-1, -1], [-1, 1], [1, -1], [1, 1]]

The same for edges.

Download PS.mw

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