restart;
M:=2; k:=1;
U:=Matrix(M+1, (r,s) -> `if`((r>1) and (s<r) and (r+s)::odd, 2^(k+1)*sqrt((2*r-1)*(2*s-1)), 0));
DD:=LinearAlgebra:-DiagonalMatrix([U $ 2^k]);

To compute the limit of a Riemann sum, the command eulermac usually works.

But your sum Sum(csc(Pi*x/i)^3*sin(Pi*x)^3/i^3, i = 1 .. n) is not a Riemann sum. And its limit for n-->infinity is actually the series

Sum(csc(Pi*x/i)^3*sin(Pi*x)^3/i^3, i = 1 .. infinity);

It's easy to see that this series diverges for each x in R \ Z (for x in Z the sum is not defined).

So, you have to rethink your problem.

The error message is clear: you must have numeric values for y, because Analytic is designed for expressions of a single variable.
Try for example:
RootFinding:-Analytic(eval(F, y=3), z, re = -5 .. 0, im = -100 .. 100);

The main difference between

f := x -> expr1;
and
g := unapply(expr2, x);

is that expr2 is evaluated, but expr1 is not.

So, for example

restart;
y:=x+1;
f:=x -> y;
g:=unapply(y,x);
y:=7;
f(10);  # 7
g(10);  # 11

For the  x[...], part,  a procedure cannot have indexed names as formal parameters. That is why, unapply automatically corrects e.g.  x[1] into a x_1.  Trying to define  f := proc(x[1]) ... end proc  generates an error.

UMG:=proc(n::posint); # Unitary Matrix Generator
local RM:='LinearAlgebra:-RandomMatrix(n, generator=rand(0. .. 1.))', S, i;
S:=eval(RM+I*RM);
Matrix(LinearAlgebra:-GramSchmidt([seq(S[..,i],i=1..n)],normalized));
end:

Digits:=15:
U:=UMG(3):
fnormal(U.U^*);  #check

restart;
with(Statistics):
p:= 9/10:
L := Distribution(ProbabilityFunction = (k -> -1/ln(1-p)*p^k/k), 'Support' = 1 .. infinity):
X := RandomVariable(L):
S:=Sample(X, 10000):
Histogram(S);

You need to know to write mathematical expressions in Maple:

(1 + a * cos(theta1*x)) * cos(theta2*x) = b;

I'll assume that theta1, theta2 are not 0 (the fact that theta2>theta1 is irrelevant).

If a=0 or b=0 then Maple can solve it easily. Try e.g.

solve( (1+a*cos(theta1*x))*cos(theta2*x) = 0, x, allsolutions );

Now, if a*b <> 0, the equation can be solved symbolically only if  theta1/theta2 is rational (it can be reduced to a polynomial equation and Maple will be able to solve it). Try:

solve( (1+a*cos(theta*x))*cos(2*theta*x) = b, x, allsolutions );

The expensive ifactors should be avoided.

P11:=proc(N)
# divisible by a prime > 11
local n:=N, q;
while irem(n,2,'q')=0 do n:=q od:
while irem(n,3,'q')=0 do n:=q od:
while irem(n,5,'q')=0 do n:=q od:
while irem(n,7,'q')=0 do n:=q od:
while irem(n,11,'q')=0 do n:=q od:
evalb(n>11)
end:

Check11:=proc(a,b);
local n;
for n from a to b do
  if P11(n) then next 
  elif P11(n+1) then n:=n+1; next
  elif P11(n+2) then n:=n+2; next
  elif P11(n+3) then n:=n+3; next fi;
  return n;
od;
true
end:

CodeTools:-Usage(Check11(1000000,2000000));

The help file mentions only  D(y)(x0) = y0 for initial conditions in dsolve.

I see  eval(diff(y(t),t),t=x0)=y0 just as a bonus. So, not a bug!

Just replace Norm(foo, 2)  with  Norm(foo, 2, conjugate=false)

It's not a bug.
expr mod p   ==>  mod is applied to the rational coefficients of the expression; expr is usually a polynomial or a rational expression (i.e. ratio of two polynomials).
So, 
10*i mod 7;   #  ==> 3*i
i mod 3;        # ==> i
12/13 mod 10;   # ==> 4,  because 4*13 = 12  (modulo 10)

with(plots): with(plottools):
display(sphere([0, 0, 0], 1), 
        plot3d(1.01, theta=0..Pi/2, phi=0..Pi/2, coords=spherical, color=red, style=surface),
        axes = framed );

f:=proc(a,b,ev)
   local i;
   if ev<>true then return "Waiting" fi;
   add(evalf(sin(a+b+i)), i=1..10^5)
end:
Explore(f(a,b,ev), parameters=[[a=0.0 .. 1.0],[b=2.0 .. 5.0],[ev, controller=checkbox]]);

Linearization (near x0) simply means replacing f(x) with f(x0) + f'(x0) (x - x0).

Note that in finite dimensions, f'(x0) is identified with the Jacobian.

So, of course Maple is able to do it, and it's simple, as you see.