The first bug is that Mode returns the set of critical points (actually float approximations) instead of the maximum point.
Trying to isolate the error, a bug in solve related to LambertW function appeared. Here is the commented Maple code.

with(Statistics):
X := RandomVariable(Normal(-2,1)):
Y := RandomVariable(Normal(2,1)):
r    := 4/10:      # r=1/2 is similar
f:= (1-r)*PDF(X,t)+r*PDF(Y,t);

f__Z := unapply(f, t);
Z    := Distribution(PDF=f__Z):
Mode(Z);
   Warning, solutions may have been lost 

maximize(f,location);  #does not work, returns unevaluated!
maximize(f ,t=-10..10,location);  # wrong answer,  the max is not at t= -2

f1:=diff(f,t):
solve(f1,t);   # WRONG!!!  returns the result of Mode
s:=allvalues(%);
eval(f1, t=s[1]);  evalf(%);  #  so, not a root!
     

        -0.2141283611e-3

ans:=fsolve(f1,t=-3..-1);
      ans := -1.999102417
evalf(eval(f1, t=%));
     

## So, the true Mode is ans  [as float approx].

In Maple only the first question has to be answered, the rest are automatic (almost).
See:
?sum   ?exp   ?infinity

Hint: start with
assume(delta <= 0);
(for convergence).

Q:=Quantile(X, z);
simplifies to
RootOf(_Z^5-30*_Z^3-875*_Z+10000*z-5000)

But this is not quite correct; the polynomial in _Z may have more than one real root and the correct one is in -5..5). The plot in the first case is for another branch. In the second case, Q is not evaluated and the plot is correct but very slow (because probably fsolve is called)!

The fast approach is to choose the correct branch in RootOf.
RootOf(op(Q), -5..5);
plot(%, z=0..1);
 

mtaylor(sin(x+y)+y, [x=Pi, y=0], 6);

sinc:= x -> sin(x)/x:
sinc(0):=1: sinc(0.):=1.:
n:=10;  #1000 if you insist
a:=5.0; h:=a/n;
Z:=Matrix(2*n+1, (i,j) -> sinc((i-1-n)*h)^2*sinc((j-1-n)*h)^2, datatype=float[8]);

# note that sinc is even, so Z is symmetric.

You can convert using IrfanView (free)
http://www.irfanview.com/main_formats.htm

Explore(plot([(-y0-x0)*exp(-2*t)+(y0+2*x0)*exp(-t),
               -2*(-y0-x0)*exp(-2*t)-(y0+2*x0)*exp(-t), t=0..5]),
        x0=1..5., y0=1..5.);

evalb(irem(390,7)=0);
   false

Q:=[X,Y]:
Ps:=proc(f,g)
local k;
global P,Q; # not really necessary
add(diff(f,P[k])*diff(g,Q[k]) - diff(f,Q[k])*diff(g,P[k]),k=1..nops(Q))
end: 

Ps(P[1],X) ;

with(plots):with(plottools):
X:=cos(t):Y:=sin(t):
vol:=Pi*X^2*(2*Y) + 2*Pi*X^2*(1-Y)/3:
sol:=maximize(vol,t=0..Pi/2,location):
t:=eval(t, op([1,1],sol[2]));#evalf(%)*180/Pi;
p1 := display(polygon([[[-X, Y], [0, 1], [X, Y], [X, -Y], [0, -1], [-X, -Y], [-X, Y]]]), color = red, thickness = 2):
p2 := implicitplot(x^2+y^2 = 1, x = -1 .. 1, y = -1 .. 1):
display(p1, p2);
simplify([X,Y]);

#[X,Y]

It is not fsolve, it is your system which contains

[c0 = 0, ... , (...)/c0 + ... = 0, ...]

Of course it cannot be solved. Replace ":" by ";" to see what happens.

PS. Why do you use such strange constructs?
F(eta) := ...
...
F(eta) := unapply( F(eta), eta);
Sintactically it is correct but now F is an empty procedure and in its remember table F(eta) is also a procedure. It works, but why so?
 

with(Statistics):
m:=Mean(a):
v:=Variance(a):
b:=v/m: c:=m/b:
P:=DensityPlot(RandomVariable(Gamma(b,c))):
plots[display](h,P);

You may clip the plots manually instead of using view.
For example,

g:=sin(k*x):
g1:=piecewise(g<1/2, g, undefined):
plots:-animate(plot, [g1, x=0..4, view=[-6..6,-2..2]], k=1 .. 6.0);

You cannot expect acceptable results for this model: your function takes values < 1 (for any parameters) while your y-vector has almost all values > 20.

In Windows you could write a short Autoit script to "export as .mpl" (see https://www.autoitscript.com) and insert a ssystem call  in the worksheet which executes the script.