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.

The result above is incorrect, since the straight line and the parabola in this example intersect not at one, but at two points. To avoid errors for different curves, use different parameter names. Here is the right solution:

Download parametrics.mw

Red arrows (tangent vectors) show the direction of increasing parameter:
 

Download arrows1.mw

Try the following:

plot([x+deltax, z, x=0..l]);
 

Here is the plot that OP wished according to his link:

restart;
f:=x->sqrt(1-x^2):
g:=x->-sqrt(1-x^2):
y1:=f(a)+D(f)(a)*(x-a):
y2:=g(a)+D(g)(a)*(x-a):
y3:=[-1,t,t=-1.5..1.5]:
y4:=[1,t,t=-1.5..1.5]:
Opt:=x=-2..2,color="Blue":
plots:-display(plot([seq(op([y1,y2]),a=-0.99..0.99,0.02)],Opt),plot([y3,y4],Opt), view=[-2..2,-1.5..1.5], scaling=constrained, size=[700,500],tickmarks=[spacing(0.5),spacing(0.5)],  axis[1]=[gridlines=40], axis[2]=[gridlines=30]);

Output:
                    

L := [8,7,9,12];
L1:=[seq([i, L[i]], i=1..nops(L))];
map(t->t[1], sort(L1, (x,y)->y[2]<x[2]));
 

Your second result obtained with erroneous syntax (numeric option is  misplaced) is incorrect. To test this, I did both calculations for  T=1, after making some simplifications to your commands (Pi should be instead of pi, and also removed some commands and options that do not affect the result). I don't know how Maple could get it but the second result is also erroneous because the integral in which the integrand is positive must also be positive.

Download int_(2)_new.mw

What function do you call monotonically increasing (or monotonically decreasing) on some multidimensional set? The issue is that there is a natural order on the number axis: for any two different points, we know for sure which one is greater. And in multidimensional space? For example, on 2D which point is greater: (0,1) or (1,0)?

I do not understand why such a complicated code (nested loops and etc.) for a simple task. You can immediately build all the necessary plots, if you know their equations, and then get the data using  plottools:-getdata  command:

restart;  
H:= 3+2*tanh(x-6);
P:=plot([H, -H, 0], x=0..10, linestyle=[1,1,3], color=black, thickness=[2,2,0], view=[-1..12,-6..6], axes=box);
interface(rtablesize=1000);
map(t->t[3],[plottools:-getdata(P)]);

restart;
f:=ln((1-x)^2*(x+1)^2/((-I*x-I+sqrt(-x^2+1))^2*(I*x+I+sqrt(-x^2+1))^2));
applyop(t->numer(t)/expand(denom(t)), 1, f);
simplify(%);
applyrule(ln((-1+x)^2)=2*ln(abs(x-1)), %);
sort(%);
                

Edit.