Symbolic pdsolve produces particular solutions, but does not handle initial and boundary conditions.  Numerical pdsolve handles initial conditions and boundary conditions, but will have trouble with yours because of the singularity at r=0.

The classic way to solve your problem is this.

Find a constant solution of the pde with boundary conditions: it turns out to be

-Beta*a/sinh((lambda/Di)^(1/2)*a)/lambda/r*sinh((lambda/Di)^(1/2)*r)+Beta/lambda

Do separation of variables on the homogeneous pde, then use the boundary conditions to get a basis of solutions: they turn out to be

exp(-Pi^2*Di*n^2*t/a^2)*sin(n*Pi*r/a)/r*exp(-lambda*t)

So you want a solution of the form

-Beta*a/sinh((lambda/Di)^(1/2)*a)/lambda/r*sinh((lambda/Di)^(1/2)*r)+Beta/lambda+Sum(c[n]*exp(-Pi^2*Di*n^2*t/a^2)*sin(n*Pi*r/a)/r*exp(-lambda*t),n = 1 .. infinity)

For t = 0 the initial condition says

-Beta*a/sinh((lambda/Di)^(1/2)*a)/lambda/r*sinh((lambda/Di)^(1/2)*r)+Beta/lambda+Sum(c[n]*sin(n*Pi*r/a)/r,n = 1 .. infinity) = 0

Multiply by r and you have a Fourier sine series.  To get the coefficient c[n], multiply by sin(n*Pi*r) and integrate r from 0 to a.  I get

> c[n] := -2/a*int((-Beta*a/sinh((lambda/Di)^(1/2)*a)/lambda*sinh((lambda/Di)^(1/2)*r)+Beta/lambda*r)*sin(n*Pi*r/a),r=0..a) assuming n::posint;

c[n] := 2*Beta*a^3*(-1)^n/(lambda*a^2+Di*n^2*Pi^2)/n/Pi

 So the solution is

-Beta*a/sinh((lambda/Di)^(1/2)*a)/lambda/r*sinh((lambda/Di)^(1/2)*r)+Beta/lambda+Sum(2*Beta*a^3*(-1)^n/(lambda*a^2+Di*n^2*Pi^2)/n/Pi*exp(-Pi^2*Di*n^2*t/a^2)*sin(n*Pi*r/a)/r*exp(-lambda*t),n = 1 .. infinity)

max should work directly on an Array (in recent versions of Maple).  No need to flatten.

max(Img_pigeon_HWT);

This is not usually called a "probability distribution function", it is a "probability mass function" or "probability function".  Moreover you want this for integers x from 0 to n.

f:= (n,x) -> n!/(x!*(n-x)!*(2^(n)));
plots[pointplot]([seq([x,f(100,x)],x=0..100)]);

Or perhaps you'd prefer

with(Statistics): 
X:= RandomVariable(Binomial(100,1/2));
DensityPlot(X);

To answer your second question, what you'd want would be the cumulative distribution function
 

plot(CDF(X,x), x=40 .. 60);

For 100 tosses, you can just add the probabilities of x heads for x from 45 to 55.

Or you could use the cumulative distribution function...

In n tosses, the probability of getting exactly x heads is binomial(n,x)/2^n, where x is any integer from 0 to n.

This is called reversion of a series. 

For example:

> AS:= add(a[n]*R^n, n=1..9);
  Sol := series(RootOf(AS - A,R), A, 10);

Sol := A/a[1]-a[2]*A^2/a[1]^3-(-2*a[2]^2+a[3]*a[1])*A^3/a[1]^5-(5*a[2]^3-5*a[3]*a[2]*a[1]+a[4]*a[1]^2)*A^4/a[1]^7-(-14*a[2]^4+21*a[2]^2*a[3]*a[1]-6*a[4]*a[2]*a[1]^2-3*a[3]^2*a[1]^2+a[5]*a[1]^3)*A^5/a[1]^9-(42*a[2]^5-84*a[3]*a[2]^3*a[1]+28*a[4]*a[2]^2*a[1]^2+28*a[2]*a[3]^2*a[1]^2-7*a[5]*a[2]*a[1]^3-7*a[4]*a[3]*a[1]^3+a[6]*a[1]^4)*A^6/a[1]^11-(-132*a[2]^6+330*a[3]*a[2]^4*a[1]-120*a[4]*a[2]^3*a[1]^2-180*a[3]^2*a[2]^2*a[1]^2+36*a[5]*a[2]^2*a[1]^3+72*a[4]*a[2]*a[3]*a[1]^3-8*a[6]*a[2]*a[1]^4+12*a[3]^3*a[1]^3-8*a[5]*a[3]*a[1]^4-4*a[4]^2*a[1]^4+a[7]*a[1]^5)*A^7/a[1]^13-(a[8]*a[1]^6+90*a[5]*a[2]*a[3]*a[1]^4+429*a[2]^7+990*a[3]^2*a[2]^3*a[1]^2-165*a[3]^3*a[2]*a[1]^3+495*a[4]*a[2]^4*a[1]^2-165*a[5]*a[2]^3*a[1]^3+45*a[4]*a[3]^2*a[1]^4-9*a[7]*a[2]*a[1]^5+45*a[6]*a[2]^2*a[1]^4-9*a[6]*a[3]*a[1]^5+45*a[2]*a[4]^2*a[1]^4-9*a[5]*a[4]*a[1]^5-1287*a[3]*a[2]^5*a[1]-495*a[4]*a[2]^2*a[3]*a[1]^3)*A^8/a[1]^15-(110*a[5]*a[2]*a[4]*a[1]^5+110*a[6]*a[2]*a[3]*a[1]^5-5*a[5]^2*a[1]^6+a[9]*a[1]^7-55*a[3]^4*a[1]^4-5005*a[3]^2*a[2]^4*a[1]^2+5005*a[3]*a[2]^6*a[1]-10*a[8]*a[2]*a[1]^6+55*a[7]*a[2]^2*a[1]^5-220*a[6]*a[2]^3*a[1]^4+715*a[5]*a[2]^4*a[1]^3+1430*a[3]^3*a[2]^2*a[1]^3-2002*a[4]*a[2]^5*a[1]^2+55*a[5]*a[3]^2*a[1]^5-10*a[7]*a[3]*a[1]^6+55*a[3]*a[4]^2*a[1]^5-330*a[4]^2*a[2]^2*a[1]^4-10*a[6]*a[4]*a[1]^6-1430*a[2]^8+2860*a[4]*a[2]^3*a[3]*a[1]^3-660*a[4]*a[2]*a[3]^2*a[1]^4-660*a[5]*a[2]^2*a[3]*a[1]^4)*A^9/a[1]^17-(4862*a[2]^9+1001*a[3]^4*a[2]*a[1]^4+2002*a[4]^2*a[2]^3*a[1]^4-11*a[6]*a[5]*a[1]^7+66*a[2]*a[5]^2*a[1]^6-3003*a[5]*a[2]^5*a[1]^3+66*a[8]*a[2]^2*a[1]^6+1001*a[6]*a[2]^4*a[1]^4-286*a[7]*a[2]^3*a[1]^5-11*a[9]*a[2]*a[1]^7+8008*a[4]*a[2]^6*a[1]^2-286*a[4]*a[3]^3*a[1]^5-11*a[8]*a[3]*a[1]^7+66*a[6]*a[3]^2*a[1]^6-19448*a[3]*a[2]^7*a[1]+24024*a[3]^2*a[2]^5*a[1]^2-10010*a[3]^3*a[2]^3*a[1]^3+22*a[4]^3*a[1]^6+a[10]*a[1]^8-15015*a[4]*a[3]*a[2]^4*a[1]^3+6006*a[4]*a[3]^2*a[2]^2*a[1]^4-858*a[6]*a[2]^2*a[3]*a[1]^5+132*a[7]*a[2]*a[3]*a[1]^6+4004*a[5]*a[2]^3*a[3]*a[1]^4-858*a[5]*a[2]*a[3]^2*a[1]^5+132*a[5]*a[3]*a[4]*a[1]^6+132*a[6]*a[2]*a[4]*a[1]^6-858*a[5]*a[2]^2*a[4]*a[1]^5-858*a[4]^2*a[2]*a[3]*a[1]^5-11*a[7]*a[4]*a[1]^7)*A^10/a[1]^19+O(A^11)

n = 30 took a few minutes.  Of course the result is too big to reproduce here: the length is 1408790.

When you say you have 60 Gb, I suspect that's disk space, not RAM, unless your computer is _really_ high-end.  While you could use that space as virtual memory, when Maple starts using virtual memory things typically slow down to less than a crawl.

All frames in an animation use the same scaling of x,y,z axes.  To zoom in, you'll want to rescale the input to mandelbrot.  Thus:

> for i from 0 to 50 do
     man:= (x,y) -> mandelbrot(-2 + .0398*i + (2.5 - .04988*i)*x,  -1 + .036078*i + (2 - .039898*i)*y);
     frame[i]:= plot3d(0, 0..1, 0..1, color=man, orientation=[-90,0],style=patchnogrid,grid=[100,100]);
   end do:
   plots[display]([seq(frame[i],i=0..50)],insequence=true);

www.math.ubc.ca/~israel/manzoom.gif

You're going to integrate over the region of the xy plane inside both of the circles (x-2)^2 + y^2 = 1 and x^2 + y^2 = 4.  That looks like this:

> X:= r*cos(t): Y:= r*sin(t): with(plots): with(plottools):
display([circle([2,0],1),circle([0,0],2),
implicitplot((X-2)^2+Y^2<=1,r=0..2,t=0..2*Pi,coords=polar,filledregions=true)],
scaling=constrained);

Let's put the x integral on the inside and the y integral on the outside. 

The endpoints for the y integral will be the y coordinates of the two points where the circles intersect.

For the x integral, the endpoints will be the x coordinates of a point on the small circle and a point on the large circle...

Besides writedata, which would produce a text file that could be imported into Excel, you could produce an actual Excel worksheet with Export in the ExcelTools package.  For example:

> A:= Matrix(5,5,(i,j) -> 1/(i+j));

> ExcelTools[Export](A, "c:/mypath/myfile.xls");

> plot3d(exp(x)+y, x = 0 .. 5, y = 0 .. 5, axis[3] = [mode = log], axes = box);

Hint: if A and B are points represented as vectors, the midpoint of AB is (A+B)/2.

What have you done so far?

Note that the region in question has its top and bottom on the sphere, and its projection on the xy plane is the intersection of two circles.
Express the volume using an iterated integral.

There is probably no closed-form solution.  All dsolve can do is reduce it to a first-order nonlinear equation.

Numerical solutions, series solutions, and phase plane analysis are possible, depending on what the rest of the assignment is.

testeq does not "show they are same in real domain".  In fact, testeq does not know anything about real numbers, and is quite unreliable for expressions involving complex exponentials. In this case g1 and q[1] are definitely not the same. 

> simplify(eval(g1-q[1],{rho=2,delta=4,R=3,a[t]=0,a[t+2]=0,inc[t]=0,inc[t+1]=1,_Z3=0,_Z9=1}));

                         6

Here is a simpler example of the trouble with testeq:

> testeq(exp(2*Pi*I*x)=1);

      true