The numerical integration routine for

is fooled by the fact that the integrand is 0 in the intervals [2k, 2k+1)  and stops the summation too early.

map(min, M, 20);

or maybe

map[inplace](min, M, 20);

You have used distinct strange names in the two files:

`#mover(mi("Ξ"),mo("&uminus0;"))`
and
`#mover(mi("Ξ",fontstyle = "normal"),mo("&uminus0;"))`;

I'd suggest to avoid in the future such names in files (you cannot see them properly in 2d input!).
The fix is simple:
 

restart;
read "sai.m";
`#mover(mi("Ξ"),mo("&uminus0;"))`:
`#mover(mi("Ξ",fontstyle = "normal"),mo("&uminus0;"))`:=%:
read "TKtm.m";
TKtm;

      0.106923339809255e-6*(diff(tau[1](t), t))^2

printlevel := 2; 
Equation := 4;
for i from 0 to Equation do 
for j from 0 to Equation do 
  C[i,j]:=coeff(coeff(H1,x,i),y,j)
end do end do;

The documentation says (see ?_):

Any symbol beginning with an underscore (_) is effectively reserved for use only by library code. It is not available to users. Failure to observe this rule can lead to unexpected results.

p:=2*yd[0]*k[a1]*k[d1]*ya[1]+(alpha*C[T]*k[a1]*k[m]-alpha*R[b]*k[a1]*k[d1]-alpha*R[m]*k[a1]*k[d1]-alpha*k[d1]*k[m])*ya[1]-2*k[a1]*k[d1]*yd[1]*yd[0]+(-alpha*C[T]*k[a1]*k[m]+alpha*R[b]*k[a1]*k[d1]+alpha*R[m]*k[a1]*k[d1]+alpha*k[d1]*k[m])*yd[1]:

sort([coeffs(p, indets(p, {ya[integer],yd[integer]}))], length)[1];

P:=[
y[a0]-y[d0],
k[d1]*y[a1]-k[d2]*y[d2],
k[d1]*y[a1]*x[1]-k[d2]*y[d2]*x[2]
]:

remove(has, subsindets(P, specindex(integer,x), _XX), _XX);

You should be aware that

plot(f(x), x=a..b, coords = logarithmic);

is equivalent to

plot( [ln(f(x)^2+x^2)/Pi, 2*arctan(x/f(x))/Pi, x=a..b]);

So, use coords = ...  only if you know exactly what it means.

[subsindets(r, {set,list}, op)];

        [a, b, c, d, e, f, g, h]

The integral is most probably divergent.

Note that
int(exp(-abs(x1-x2)), [x1 = -infinity .. infinity, x2 = -infinity .. infinity]) = infinity

It seems that Mathematica knows about the recurrence  f(n+1) = M(f(n)), where M is a Moebius function, but Maple does not.
(It should be easy to implement it.)
Workaround:

A:=LinearAlgebra:-MatrixFunction(<1,1;2,0>, x^(n-1), x):
V:=A.<-7/13,1>:
simplify(V[1]/V[2]);

PDF has severe limitations. For example the PDF for a sum of two random variables fails if the random variables are not independent.
A more concrete example: if X,Y are random variables obtained via RandomVariable(Uniform(0,1))
then I have noticed that

PDF(X^2 + Y^2 + a*X*Y, x);

works only for a ∈ {0, 2, -2}.

# f(x,y) = h(x*y), where
# h := z -> minimize(sin(a*z), a=1..5);
# But h(z) > 0 iff z<Pi/5,  so, your int equals
int(Pi/(5*x), x=Pi/5..1) + Pi/5; evalf(%);

       (1/5)*Pi*ln(5)-(1/5)*Pi*ln(Pi)+(1/5)*Pi
       .9203031946

The 3rd line actually does not solve the ode. It just reduces it to a second order ode. Practically it only makes the substitution y'' = z.

The NULL answer for the 4th line is normal. dsolve cannot solve it for a general K(x); the solution may very well not exist for many K(x), P and L.