You have a hard to detect typo:

ics := Y(b) = Ts-Tb;  (D(Y))(-b) = 0;

instead of

ics := Y(b) = Ts-Tb,  (D(Y))(-b) = 0;

In the first place you should not unprotect GAMMA (which is an important function in Maple). Why not use Gamma instead which is free and looks the same?

You also use ArrayTools:-Alias which is used for low-level manipulation and could be incompatible with Grid (I have not tested).
See also the help page of Grid:-Seq about Grid:-Set.
Do you obtain the expected results when using seq instead of Seq?

P.S. Using 2D math input + Document mode + "special" names for variables such as your `#mover(mi("Γ",fontstyle = "normal"),mo("&uminus0;"))`   makes debugging very difficult.  This could be OK for a presentation but not for programming!

You want the singular points.
For your
expr := (a-b) / (b-a^3);    # be careful with the parentheses!
singular( expr, {a,b} );
      {a = a, b = a^3}

So, you have infinitely many singular points (a is arbitrary and b=a^3).

It is not clear what are your intentions. For a minimization problem, usually the command detects the singularities.
E.g.

minimize(1/x, x=-1..2, location);
    -infinity, {[{x = 0}, -infinity]}

You cannot. And you do not need this.
CompressedSparseForm  is supposed to be used for external programs/platforms.

L must be numeric.

restart;
with(IntegrationTools):
L:=5:
simplify(int(f(x), x = 0 .. L*Ts)-Split(int(f(x), x = 0 .. L*Ts), [i*Ts, i = 0 .. L])):lprint(%);
       int(f(x), x = 0 .. 5*Ts)-(int(f(x), x = 0 .. Ts)) -
         (int(f(x), x = Ts .. 2*Ts))-(int(f(x), x = 2*Ts .. 3*Ts))-(int(f(x), x = 3*Ts .. 4*Ts))-(int(f(x), x = 4*Ts .. 5*Ts))
combine(%);

         0

If you just want to parallelize some numerical computations (without int, sum etc, see ?thread-safe) you can use Threads:-Seq.
Example

f := proc(i) local k; 
if   i=1 then    add(1/evalf((k+sin(k))^(k/(k+1))),k=1..5*10^4)
elif i=2 then    add(1/evalf((k+sin(k))^(k/(k+2))),k=1..6*10^4)
elif i=3 then    add(1/evalf((k+sin(k))^(k/(k+3))),k=1..6*10^4)
elif i=4 then    add(1/evalf((k+sin(k))^(k/(k+4))),k=1..3*10^4)
fi
end:

CodeTools:-Usage([Threads:-Seq(f(i), i=1..4)]);
CodeTools:-Usage([seq(f(i), i=1..4)]);  # without threads, for comparison

P:=proc(n::nonnegint,k::nonnegint) 
if n<1 or k<1 or k>n then return NULL fi;
if k=1 then return [n] fi;
if k=n then return [1$n] fi;
seq([u[]+~1], u=[P(n-k,k)]), seq([1,u[]], u=[P(n-1,k-1)])
end: 

CodeTools:-Usage(nops([P(100, 5)]));
memory used=323.13MiB, alloc change=76.01MiB, cpu time=2.07s, real time=1.94s, gc time=374.40ms
                             38225

CodeTools:-Usage(nops([P(80, 10)]));
memory used=3.28GiB, alloc change=435.74MiB, cpu time=29.73s, real time=23.37s, gc time=9.70s
                             533975

sort~(MM);          

1. The exponential of a vector v does not exist in maths. If you want to apply exp elementwise then use exp~(v)  or map(exp, v)
==>  < exp(v[1]), exp(v[2]), ... > .  Note that for a square matrix M there is MatrixExponential(M).

2. To integrate a matrix M (elementwise also)  use  map(int, M, t=a..b)  or map(int, M, [x=a..b, y=c..d,...] ).

3. You cannot have in an integral both  vx  and  vx(x,t)  and integrate wrt  vx.

LinearMultivariateSystem has a bug.
Use solve or SolveTools:-SemiAlgebraic  instead.

P.S. If you really need LinearMultivariateSystem  then:

sys1:=subsindets(sys, indexed, e -> cat(op(0,e),__, op(e)));
SolveTools:-Inequality:-LinearMultivariateSystem(sys1,indets(sys1));

This form of eval is an "intelligent" subs. In the second example, a subs was possible. In the first one, subs would produce a nonsense, so the expression was left practically untouched (you can check this using lprint). It is the typesetting system which produced the ...|{x=0} display.

NonegComb:=proc(M::Matrix,v::Vector)
local n:=upperbound(M,2), L, t, S;
if upperbound(M,1) <> upperbound(v) then error "Incompatible vector" fi;
L:=convert(add(M[..,i]*t[i],i=1..n)-v,list);
S:=simplex:-minimize(0, [(L=~0)[],  seq(t[i]>=0, i=1..n)]);
if S={} then return false fi;
eval([seq(t[i],i=1..n)],  S);
end:

M:=Matrix([[3, 0, 2, 2, 1, 1], [-1, -1, 0, -1, 0, -1], [0, 3, 0, 1, 1, 2], [3, 0, 1, 2, 0, 1], [-1, -1, -1, -1, -1, -1]]):
v1:= Vector[column](5, [3, 1, 0, 0, -2]):
v2 := M[..,1]+7*M[..,2]+9*M[..,4]:
NonegComb(M,v1);
                             false
NonegComb(M,v2);
                     [0, 13/2, 0, 21/2, 0, 0]

P:=n->seq(seq([a,b,n-a-b],b=max(a, floor(n/2)-a+1)..floor((n-a)/2)), a=1..floor(n/3)):

P(36); nops([%]);
[2,17,17],[3,16,17],[4,15,17],[4,16,16],[5,14,17],[5,15,16],[6,13,17],[6,14,16],[6,15,15],[7,12,17],[7,13,16],[7,14,15],[8,11,17],[8,12,16],[8,13,15],[8,14,14],[9,10,17],[9,11,16],[9,12,15],[9,13,14],[10,10,16],[10,11,15],[10,12,14],[10,13,13],[11,11,14],[11,12,13],[12,12,12]
      27