I can't understand what you want, and you can't understand what I've given you. Please get someone to help you with your English before posting.

You wrote:

The Maxwell-Boltzmann distribution sampling times are 2002. But you just list five sequences. I don't know why.

Five was just a convenient number that fit on one line on the screen. The way that you are using the word "sampling" is not standard mathematical English. The entire population of possible distributions of 9 quanta among 6 particles has size 2002. But that is neither a sample nor a time.

I am very confuced about the sampling. How can you set sampling times are 2002?

I have no idea why you would want a sample of size 2002. The idea seems crazy to me. But you can change the five to any number that you want. That should be obvious!

'R()' $ 2002;

Now, if you want the whole population, that's another matter. Now (in the worksheet posted below) I have included an option to return the whole population: ouput= population.

You said the computation is exact---not based on random sampling. But latter you said that does either the exact computation or provides a random sampler depending on an input option.

If you say output= exact then the computation is exact! If you say output= sampler then the output is a random sample! Why can't you understand that? Please get someone to help you read the code. The computation where I compared the numbers with those on the website were the exact computations.

What's the difference from random sampling and random sampler?

I am using "sampler" to mean a procedure that generates a random sample. When you call my procedure with the option output= sampler then the object returned is itself a procedure---a procedure which draws one item from the population each time it is called.

Download QuantumDistribution.mw

@H-R It is usually easier to specify a unique characterization of the component that you want to extract than it is to fix its operand number. For your situation, could I characterize the part that you want to extract as "the factor in the sum that has a sinh in it"? Would that be general enough to cover every situation where you would use this and yet precise enough to uniquely specify the part that you want?

For the case at hand, the following extracts the subexpression that you want:

select(hastype, op(1, indets(BC4, specfunc(anything, sum))[]), sinh(algebraic));

This assumes that BC4 has only one sum.

m:= 1:  N:= 2:
Z:= Matrix(m+1, m+1):
Id:= Matrix([1$m+1], scan= band[0]):
A:= Matrix([[Z$N], [Id$N-1]], scan= band[0,1]);

See ?scan (unfortunately there are no examples on this help page).

Here's another solution completely in polar coordinates. Unlike Kitonum's solution (which I don't mean to disparage), I derive the polar form of the necessary formulae for slope and tangent line from first principles.

restart:
R:= theta-> 3 - 3*cos(theta):
theta0:= 3*Pi/4:
polar:= (r,theta)-> {x = r(theta)*cos(theta), y(x) = r(theta)*sin(theta)}:
dy_dx:= PDEtools:-dchange(polar(r,theta), diff(y(x),x)):
line:= PDEtools:-dchange(polar(r,theta), y(x) - eval(y(x), polar(R,theta0)) =
     m*(x - eval(x, polar(R,theta0)))):
line_polar:= solve(line, r(theta)):
m:= eval(dy_dx, [r= R, theta= theta0]):
P1:= plots:-polarplot(R(theta), theta= 0..2*Pi):
P2:= plots:-polarplot(line_polar, theta= Pi/4..Pi):
plots:-display([P1,P2]);

Here's a procedure for it:

Pbar:= proc(P::Matrix(square), k::posint)
local
     m:= op([1,1], P),
     km:= k*m,
     M:= Matrix(km, km, (i,j)-> `if`(i < j, 1/m, 0)),
     n
;
     for n to k do M[(n-1)*m+1..n*m, (n-1)*m+1..m*n]:= P end do;
     M/k
end proc:

Example of use:

Pbar(< 1,2 ; 3,4 >, 3);

One way is by declaring the parameter as type evaln:

mysum:= (data::evaln(list))-> printf("The sum of %a is %3.3f.", data, `+`(eval(data)[])):

listA:= [1,2,3]:  listB:= [4,5,6]:
mysum(listA);
The sum of listA is 6.000.

Note that the parameter itself refers to the assigned name and that the contents of the name are accessed by using eval.

You are missing a left square bracket before b1, and Gamma should be GAMMA.

If that doesn't solve your problem, then I need to know what version of Maple you are using.

When you receive a "lost connection to kernel" error, you need to

1. save your worksheet
2. close your worksheet
3. open it again

Okay, here is one procedure that does either the exact computation or provides a random sampler depending on an input option.

THe result for Bose-Einstein energy level 4 given on that website must be wrong. If you multiply the energy levels times their average values and add them up, you must get 9, right? If you do that for the numbers in the Bose-Einstein column on that website, you will see that you don't get 9. My results come from exact computation on the entire population of possible states.

Download QuantumDistribution.mw

The symbol has no mathematical meaning in Maple, but for typing in mathematical text mode, the ordinary vertical bar works fine |. On my keyboard, it's the shifted version of the backslash key.

P1:= PlanePlot(-3*x+ 2*y + 3*z = 0, [x,y,z], normaloptions= [shape= harpoon], showbasis):
P2:= PlanePlot(
     2*y+3*z = 0, [x,y,z], normaloptions= [shape= harpoon], showbasis,
     planeoptions= [color= red]
):
plots:-display([P1,P2], caption= "");

Here are two ways to accomplish the task with rtable_scanblock. This first one is in the spirit of the examples given at ?rtable_scanblock :

ind:= rtable_scanblock(A, [], (val,ind,res)-> `if`(val >= 3, [ind,res[]], res), []):

That's pretty good timewise and has the benefit of being a one-liner. Timewise it's about midway between Preben's Matrix solution and Pagan's seq solution. But it suffers a bit from the res being constantly appended and passed back and forth. The following is 3 to 4 times faster on the 2000x2000 Matrix that I used for testing. Instead of using res for the results, I store the results in a table.

Ind:= table():
rtable_scanblock(A, [], (val,ind)-> if val >= 3 then :-Ind[ind]:= 0 fi):
ind:= [indices(Ind, nolist)]:

Note that the value 0 that I assign is insignificant; it could be any value. All that matters are the indices of the table that are assigned.

Update: The following variation of Joe Riel's dynamic Array solution (given in a subsequent Answer) fares about the same timewise as my table-based solution:

Ind:= Vector():  k:= 0:
rtable_scanblock(A, [], (val,ind)-> if val >= 3 then :-k:= k+1; :-Ind(k):= ind fi):

The commands diff and convert are not thread safe. See ?index,threadsafe .

The following procedure computes all three distributions. The computation is exact---not based on random sampling. The numbers match perfectly with those on the web site that you linked to.

Download QuantumDistribution.mw

Although I strongly prefer Markiyan's use of display(..., insequence= true) to my use of delay-evaluation quotes below, I wanted to point out what was wrong with the original animate command since the problem is quite subtle. Once the call to P(n) is made, defining the parameterized random variable, it is impossible to substitute a numeric value for n. This is a limitation, severe in my opinion, of the Statistics package. But you can work around it by using delayed evaluation in the animate command like this:

animate(DensityPlot, ['P(n)', range= 0..10], n= 1..5, frames= 5);

The reason that P(n) does not cause a problem in Markiyan's version is that the special evaluation rules of seq act like delayed evaluation.

On the line where you define o[1,6], you have an angle of 270. That should be 3*Pi/2.

The mathematical constant Pi is spelled in Maple with an uppercase P and lowercase i. A variable spelled pi is just another variable.

Do not use a decimal when you know the exact value of a constant. In other words, use 1/2 instead of 0.5.

After the line defining Ldq:= Tmatrix.L6.Tmatrixinv, put the following:

convert(%, rational): #Correct decimals
subs(pi= Pi, %): #Correct pi
map(expand, %):
map(combine, %);

The resulting 3x3 matrix takes 6 lines on my screen. Let me know if you think that is simplified enough.