For d=0, f is a PDF  (actually a Negative binomial distribution).

It seems to me (after some investigations) that this is the only case when f is a PDF.
According to Maple, for [r=3, b=2, d=0.01411053479]  f seemed to be a distribution
just because  Sum_x f(x)   is very close to 1 for a small d.

Why don't you use 1D math? Later, you may convert to 2D if you want. All these problems will disappear.

I think that 1D for input and 2D for output is the best choice for now. 

se:=proc(f::procedure,a::realcons,b::realcons)
local x,s;
s:=(f(b)-f(a))/(b-a)*(x-a)+f(a);
print(plot([f(x),s], x=a..b));
s;
end proc:

f:=x -> x^2:

se(f,-1,2);

rad:=u -> `*`(op(map( z->op(1,z),factors(u)[2])));

g:=expand((x*sqrt(2)+y*sqrt(3)+z*sqrt(6)+x*y*sqrt(6)+x)^6*(x+2*y)^2):
rad(g);

If _C1 is constant just use _C1 instead of _C1(t).

If _C1(t) etc are already in an expression, use

subs( [_C1(t)=_C1, _C2(t)=_C2, _C3(t)=C3], expr);

You shoud not write in the Program Files directory.
Put your library in a directory where you have write permission. It will work.

L1toL2:=proc(L1,L2)
local f;
  f:=proc(x,L) local t; member(x,L,t); t end;
map(f,L2,L1);
end;

# it works for nondistinct elements, provided that {op(L2)} subset {op(L1)}  (of course)
L1:=[a,b,a,a]; L2:=[b,b,b,a];

{op(L2)} subset {op(L1)};
                              true
P:=L1toL2(L1,L2);
                          [2, 2, 2, 1]
L1[P],L2;
                   [b, b, b, a], [b, b, b, a]

It seems that you want to approach Riemann hypothesis using Maple. :-)

Mathematica seems to have a much better algorithm to approximate Zeta(z) for a large |z|.

Maple computes Zeta(0.6+I*1000000) very slowly (many seconds!). Other CAS-es even refuse to approximate.
It would be interesting to know the algorithm used by Mathematica and its robustness.

The probability of having 0 branch breakings after 40 steps is

1-(1-0.02)^40 = 0.5542995960 (> 1/2)

If you want a simulation to compute the number of branch breaking, you should run your code N times (in a loop) and compute the average number of breaks. (N = 10000 or something).

map(factor,A);

or

factor~(A);

Replace Flux by

J:=Flux(vv, Sphere(<0,0,0>, r), inert);

to see the integral.

Maple cannot compute it as it is. But if you change the order of integration (theta first, then phi)
(by hand or using IntegrationTools) ==>

Maybe in the near future Maple will do this automatically.

Yes, the LaTeX export is hard to edit.

I use only the latex(...) command  and paste the output in the .tex file. It also must be edited because it looks worse than the Maple output.
It is also possible redirect latex(...) to a file, see ?latex

uu:=pds:-value(u(x,t), output=procedurelist);
vv:=pds1:-value(v(x,t), output=procedurelist);
plot( [ 'log(rhs(uu(1,t)[3]))', 'log(rhs(vv(1,t)[3]))', t=0..1]);

Maple has commands for Sturm polynomials.

P :=x^6-4*x^3+x-2:

ss:=sturmseq(P, x);

sturm(ss,x,-3,3); #number of roots in (-3,3]

         2

You can inspect the source code, e.g.:

interface(verboseproc=3);
print(sturmseq);
print(`sturm/sturmrem`);
print(sturm);

It also has the command

realroot(P,10^(-4));

for isolating the real roots (based on Descartes' rule instead of Sturm).

co:=coeffs(p1-p2, [ln(E),ln(K),ln(L)]);

seq(u_=0, u_=[co]);