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.

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: