Are you interested in simulating the result ?

Here's a crude monte-carlo simulation, for 3 dimensions:

t:=time():
N:=100000:
with(Statistics):
RV := RandomVariable(Uniform(0, 1)):
s:=0:

X[1]:=Sample(RV,N):
X[2]:=Sample(RV,N):
Y[1]:=Sample(RV,N):
Y[2]:=Sample(RV,N):
Z[1]:=Sample(RV,N):
Z[2]:=Sample(RV,N):
for i from 1 to N do
    s:=s+sqrt((X[1][i]-X[2][i])^2 +(Y[1][i]-Y[2][i])^2+(Z[1][i]-Z[2][i])^2);
end do:
s/N;

time()-t;

This compares well with the exact result I posted earlier.

For higher dimensions, we can easily adapt the code for an arbitrary number 
of dimensions, but obviously it will take loner to run . So, for d-dimensional space:

RV := RandomVariable(Uniform(0, 1)):
s:=0:
d:=5:
N:=10000:
for dim from 1 to d do
    A[dim]:=Sample(RV,N);
    B[dim]:=Sample(RV,N);
end do:

for i from 1 to N do
    s:=s+sqrt(sum((A[dimm][i]-B[dimm][i])^2,dimm=1..d))
end do:
s/N;

1/105 * (4+17*sqrt(2)-6*sqrt(3)+21*log(1+sqrt(2))+84*log(1+sqrt(3))-42*log(2)-7*Pi);

                                0.6617071825

http://www.math.kth.se/~johanph/habc.pdf

http://www.math.kth.se/~johanph/Eabc.mw

f3 := (n) -> (2*n)!*(lambda)^n/(((n)!^2)*2^(2*n)*(lambda+1)^(n+0.5));
lambda:=10;
Statistics[Histogram]([seq(f3(n),n=0..20)]);

You are passing an equation, not an expression to numer. Perhaps you have parentheses in the wrong place ?

Should the penultimate line be

fsolve(numer(algsubs(alpha = xi[j+1], `finθ`[j])) = 0, beta)

or (if you need complex results)

fsolve(numer(algsubs(alpha = xi[j+1], `finθ`[j])) = 0, beta, complex)

instead of

fsolve(numer(algsubs(alpha = xi[j+1], `finθ`[j]) = 0), beta)

?

In the example from the help file, you have combined the input and output variable into one, so there is a missing variable, but maple seems to have accepted it (maybe because they are combined into the 2nd variable as a list).

But in your code, Maple thinks your linearization point is the output variables, hence the error that the 4th argument is missing.

So I think you need to correctly specify and input and output variables.

If I've understood your question correctly, I think the relation you cite applies only to a Poisson process. If events occur according to a Poisson process with parameter m, then the number of events in a period has a Poisson distribution with mean m, and the time between events has a exponential distribution with mean 1/m .  This is not true when the underlying dstribution is not Poisson.

This is an interesting question. I don't know the answer. Iit requires the inverse operation of finding the mgf from the pmf. Finding the mgf from the pmf is easy enough: Sticking with a standard Poisson variate, we have

pmf:=exp(-lambda)*(lambda^j)/factorial(j):

mpf:= sum(exp(t*j)*pmf,j=0..infinity);

                    exp(lambda (exp(t) - 1))

So how can we do this in reverse ? 

BTW, for continuous random variables, I think this requires taking the inverse (two-sided) laplace
transform of the mgf.

Derivatives wrt the spatial variable m appears only in the first equation, which is of 2nd order, so you need 2 boundary conditions for that.

You also need two initial conditons, since each equation is first order in time.. 

I also don't get the error in Maple 13 or Maple 15.

Instead I get no output, which usually means Maple found no solution so it gave up, which seems correct, as per Alec's answer above - the ICs are not consistent.

pdsolve is limited to two independent variables. Also, you may need to provide additional inputs such as a hint.

But in general, it can solve inhomogenous problems

But, for example, 

PDE := x*(diff(f(x, y), y))-(diff(f(x, y), x)) = A*cos(x);

pdsolve(PDE);
                                       /1  2    \
              f(x, y) = -A sin(x) + _F1|- x  + y|
                                       \2       /

Specifically for Poisson's equation , you might find this useful:

http://www.math.unb.ca/~seahra/resources/M4503/maple_2011/poisson_eq.pdf

Could you post your code, or your worksheet, otherwise it's hard to give specific advice.

You could try here

http://www.cs.umd.edu/~oleary/LBFGS/#maple

You need to use -y(x) not -y

Yes, pdsolve is limited to two independent variables.

Perhaps you can try an integral transform. What are your initial/boundary conditions ?

Maple thinks you have more than one Theta - Theta(r,t) and Theta(z,t). Shouldn't Theta be a function of r,z and t ? 

PDE := diff(Theta(r,z, t), r, r)+(diff(Theta(r,z, t), r))/r + diff(Theta(r,z,t),z,z) = (diff(Theta(r,z, t), t))/alpha;