@H-R 

I wanted the code to be repeatable. Thus, it is necessary to delete the file if it already exists. However, if the file does not exist (such as on the first run of the program), the Remove command will give an error. The try ... catch: ... end try command "traps" the error, so that it can be ignored. The code which may generate the error comes immediately after the try. If there is something special to do when the error occurs, that comes after the catch. In this case, there is nothing that needs to be done. However, the syntax requires the presence the word catch: anyway. See ?try .

@Stephen Forrest 

It would be nice to be able to run Part1.mw more than once, i.e., if MyNewArchive.mla already exists.

@Markiyan Hirnyk 

You didn't execute the code that defines result. That is clear from the label on the horizontal axis. This is the plot that you should get

What sum knows about this seems very basic to me. Consider this generalization of your example:

sum((-1)^k*(f(k) - a), k= 1..n);

@Markiyan Hirnyk 

Your code is remarkably simple, and I'd wager that it'd come in second place timewise among the current contenders. 

@derekcowley 

See ?simplify,siderels . The command uses a Groebner basis. The efficiency of that famous algorithm is notoriously sensitive to various factors, most notably the "monomial ordering."  It is however a true algorithm: It is guaranteed to give a result fully reduced with respect to the side relations given sufficient time and memory.

Note that your side relations are eight cubic polynomials each with nine terms. Try reducing that complexity by using left sides that are pairs of factors rather than triples. Or apply the side relations one at a time (in a series of simplify commands).

@Joe Riel 

Joe: +1 for the syntax in P do f(M) end do. It's surprising that no explicit call to g is needed in the loop. That being the case, pointers h and g are not needed. So the code can be further abbreviated to

(m,n):= (2,2):
P:= Iterator:-MixedRadixTuples([2 $ m*n]):
M:= ArrayTools:-Alias(ModuleIterator(P)[2](), [m,n]):
in P do printf("%{}d\n\n", M) end do:

Markiyan: It looks like you don't have a compiler in Maple 16. Change your call to MixedRadixTuples to MixedRadixTuples(..., compile= false): Doing this, I got it to work in Maple 16.

@Aakanksha 

It's simple:

plots[logplot]([result1, result2, result3]);

The reason that Markiyan's suggestion did not work is that your Maple is a very old version. What version is it?

@Aakanksha 

That should be plots[pointplot](result1), not plots[pointplot(result1)].

@quann169 

Please see the response by Alejandro Jakubi above and follow the link that he provided. That link leads to a post that has a video that explains how to deal with this situation. Unfortunately, you cannot use $include in Maple IDE, but apparently there is a viable alternative.

@k20057 5 The method that you describe for converting a Maxwell-Boltzmann sample to a Bose-Einstein sample does not make sense to me when you consider a large sample. For 9 quanta and 6 particles, the Maxwell-Boltzmann model has a population of 2002 different states and the Bose-Einstein model has 26 different states. If you have a sample size of 100,000, then every Maxwell-Boltzmann state is occupied about 50 times. If you convert that to a Bose-Einstein sample in the method that you describe, then all 26 states will be used with probability very close to 1. You might as well just look at the entire population of 26 states rather than a sample.

In other words, after making the conversion as you describe, there is no longer anything random about the sample.

What would make sense to me is taking the Maxwell-Boltzmann sample and reducing it by removing any members that are not in the Bose-Einstein population, and then further reducing it by removing any members that are not in the Fermi-Dirac population.

@H-R I was using Maple 18. The great increase in speed that I experienced was due to Maple 18's hardware float evaluation of Bessel functions.

Where is the N that you want to set to 100? I found an N in your code, and I set it to 100 with no problem. I also set k to 100 with no problem.

@quann169 I don't know anything about IDE, but did you remember to put the $ in the first column?

The code in your file is a statement defining a procedure. When that statement is executed by the read, it defines the procedure at the top level rather than as part of the module.

@derekcowley Oh, I think that I understand now what you want: You want to express a polynomial as a sum of terms each of which is factored a particular way. That's much more difficult than just factoring. I can't immediately think of a way to do that. Certainly the factor command cannot do it. It may be possible using simplify with side relations. The side relations would be p1 - p2 = p12, p1 - p3 = p13, etc., (every possible pair of p's is equated to a unique variable). I may be able to work on it tonight.