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.

Here is an inefficient implementation of what you want:

p:= x^3+2*x*y^2-2*y^2+x:
n:= 3:
[seq(seq(coeff(coeff(p, y, j), x, n-i), j= 0..i), i= 0..n)];

To make it more efficient, you need to use PolynomialTools:-CoefficientList, but it is somewhat of a challenge to make it conform to your very nonstandard indexing scheme. If your polynomials tend to be short, as in your example, and fairly dense with respect to your indexing scheme, then the difference in efficiency will be small.

This will only work in Maple 18+. Let me know if you need a solution for earlier Maple.

# Return the angle of a point in interval [0, 2*Pi).
Angle:= proc(P::[realcons,realcons])
local arg:= evalf(arctan(P[2],P[1]));
     `if`(arg < 0, evalf(arg+2*Pi), arg)
end proc:

SortByAngle:= (A::Matrix)-> Matrix(sort(convert(A^%T, listlist), key= Angle))^%T:

M:= < -1,0,0,-1 ; 0,1,-1,0 >:
SortByAngle(M);

The restart command should always go in its own execution group. There are various subtle things that get messed up if it isn't on its own. This is documented at ?restart . If you re-execute your code with restart in its own execution group, then the tilde will not be displayed.

Wow, it's a very interesting example. And you seem to have the same programming philosophy as me. Anyhow, you can evaluate your expression like this:

eval(subs(z= 0.5, laplace= inttrans[laplace], r));

Here it is. Please let me know if this is correct.

Am:= m-> Matrix(m+1, m+1, (i,j)-> (-1)^(j-i)*binomial(m,i-1)*binomial(m-i+1, j-i)):    
Am(9);

Please stop using that module version of RandomCompositions. Use the procedure version that I posted about an hour after I posted the module. I even explained to you in detail how the procedure worked.

RandomCompositions:= proc(n::posint, k::posint)
local
     C,
     Compositions:= [seq(C-~1, C= combinat:-composition(n+k, k))],
     Rand:= rand(1..nops(Compositions))
;
     ()-> Compositions[Rand()]
end proc:

R:= RandomCompositions(8,6):
n:= 2^13:
S:= 'R()' $ n:
T:= map(lhs=rhs/n, Statistics:-Tally(op~([S])));

plot([op]~(T), style= point, symbolsize= 16);

Your sys is a module, so start with

exports(sys);

We see that tf is an export. Each module member is accessed with the :- operator. So,

sys:-tf;

The resulting object is a Matrix, so, finally

sys:-tf[1,1];

eval([seq(y=45*r*t, t= 2..150)], r= 5);

How about this?

subsindets(map([op], [op](ex_expr1)), indexed, h-> [op(0,h), [op](h)]);

Put the legend option on the individual plots. Then the display will merge them.

restart:
with(DynamicSystems):
with(plots):
sys := TransferFunction(1/(s^2+0.2*s+1)):
p1:=ResponsePlot(sys, Step(),duration=50,color=red, legend= "step"):
p2:=ImpulseResponsePlot(sys,50, legend= "impulse");
display([p1,p2],axes=boxed, title=`step and impulse reponses`);

Suppose that you call your C with symbolic arguments, so that it returns unevaluated:

C(n,m);

If you use C:= binomial, then the output will be binomial(n,m). But if you use alias, then the output will be C(n,m). This is very useful when you use an alias for a lengthy RootOf expression---one of its common uses.