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.

eq := t = indets(f, sqrt)[];
factor(eval(f,solve(eq,{k}))) assuming t>=0;
eval(%,eq);

It's easy to setup the problem, but it seems that Maple needs a very long time to solve it.
I wonder whether it can be done in a reasonable amount of time (days?).

restart;
we:=[1,2,8,9,15,16,22,23,29,30]:
wd:=[3,4,5,6,7,10,11,12,13,14,17,18,19,20,21,24,25,26,27,28]:
cond1 := seq(seq(  add(x[i,j,k], i=1..20)>=4,  j=1..30), k=1..2):
cond2 := seq(seq(  add(x[i,j,k], j=wd)>=4,   i=1..20), k=1..2):
cond3 := seq(seq(  add(x[i,j,k], j=we)>=2,   i=1..20), k=1..2):
cond4 := seq(seq( 2*x[i,j,2]+x[i,j+1,1]+x[i,j+1,2]<=2, i=1..20), j=1..29):
cond5 := seq(p[i] = add(add(x[i,j,k],j=1..30),k=1..2)-12,  i=1..20):
cond6 := seq( [p[i]<=q[i], -p[i]<=q[i]][], i=1..20):
Optimization:-LPSolve( add(q[i],i=1..20), {seq(cond||i, i=1..6)}, binaryvariables=indets([cond1]) );

Probably you want
numapprox:-infnorm(y-mExact, 0..1);
        2.045889353*10^(-7)

Rewrite the recurrence:

rsolve({f(n+1)-f(n) = 1+f(n), f(1)=1}, f(n))

      2^n-1

Not all mathematical functions have evalhf implementation. pochhammer is one of them. See ?evalhf/fcnlist

You will have to use evalf instead.