Answering your second question, there is the command interface(prompt= ...). For example, I have the following line in my initialization file:

interface(prompt= "(**)");

Erik wrote:

I am working on a problem, which basicly is dealing with putting n identical marbles in k different boxes. I need to run through all possible combinations and do something with it.

Your "combinations" are a combinatorial structure called "compositions".

I have heard that a recursive procedure sometimes can be transformed into an iterative procedure, which will maybe make it execute faster.

That is for a small but important subset of recursions called "tail recursion". Your procedure is not that type. An example of tail recursion is

Factorial:= (n::nonnegint)-> `if`(n=0, 1, n*Factorial(n-1));

The issue with recursion is usually memory more than time.

I can't find a way to do it iteratively....

I found two pre-coded iterative ways in Maple. The first uses a very old package called combstruct. The second uses Joe Riel's (excellent) Iterator package available for free download at the Maple Applications Center here. The problem is discussed in Donald Knuth's (excellent and long-awaited) The Art of Computer Programming: Volume 4. But these methods are much slower than your recursive method. The reason to use an iterative method is that with a recursive method you essentially have to hold all combinations in memory at the same time. This is no problem for the binomial(45,5) ~ 1.3 M combinations in your example. But if you just increase your number of boxes to 10, then you have binomial(49,9) ~ 2 G, each one being an Array with 10 elements, and it's not possible in the memory of most computers.

With Joe's Iterator package, the code for your procedure is simply:

Combi:= proc(n::posint, k::posint)
local boxContent;
     for boxContent in Iterator:-BoundedComposition([n$k], n) do
          # Do something here
     end do
end proc;

With combstruct, it is

Combi:= proc(n::posint, k::posint)
local
     Iter:= combstruct[iterstructs](Composition(n+k), size= k),
     NV:= Iter[nextvalue],
     boxContent
;
     while not Iter[finished] do
          boxContent:= NV() -~ 1;
          # Do something here
     end do
end proc;

If we need to stay with this recursive procedure, then is there a way to optimize the code?

Like I said, your procedure is already quite efficient. I can make a factor-of-five improvement by using evalhf, but this will only be possible if your "# Do something here" is doable within evalhf, which is a very restricted environment. Here is your procedure done with evalhf:

Combi:= module()
export boxContent;  #A vector whose i'th entry holds the number of marbles in the i'th box

local
    Combi:= proc(k::float, boxContent::Vector, activeBox::float, marblesRemaining::float)       
    local i::integer;
        
        if activeBox < k then
            for i from 0 to marblesRemaining do
                boxContent[activeBox]:= i;
                thisproc(k, boxContent, activeBox+1, marblesRemaining-i)
            end do
        else
            boxContent[activeBox]:= marblesRemaining;
            #Do something here        
        end if;
    end proc,

    ModuleApply:= proc(n::posint, k::posint)
        boxContent:= Vector(k, datatype= float[8]);
        evalhf(Combi(k,boxContent,1,n))
    end proc
;
end module;

The call for all three procedures is simply

Combi(n, k)

I was able to do it in under a second by using simplify with side relations (see ?simplify,siderels ). I refer to it here by its syntax skeleton simplify(..., [...]). Unlike regular simplify---which may use heuristic ideas of what is simpler---simplify(..., [...]) uses a pure polynomial algebraic algorithmic approach. 

Replacing a common subexpression with a new variable is a type of side relation. Here are the substitutions that you wanted to apply:

There are a number of factors that can make a vast difference in the amount of time and memory that the algorithm uses. These factors are referred in the help page as "monomial orderings", but the help page fails to mention just how important these are. The naive approach to your problem would be

simplify(kappa, [newpar]);

I let this run several minutes until I was nearly out of memory (it was using about 4Gig) before I killed it. Then I tried applying the side relations individually from most complicated to least (this last point is very important), like this:

newkappa:= kappa:
for k from nops([newpar]) by -1 to 1 do
     newkappa:= simplify(newkappa, [newpar[k]])
end do:
newkappa;

This works nearly instantaneously (0.140 secs) and achieves a great simplification: The only variables in the result are the s's. I also did the C=4, K=7 case in 0.656 secs. So, there you go: I took it from several days of computation to under a second.

Here's the modified worksheet:

sC.mw

In the line before the else, there appears #Iota.... The # comments out the rest of the line, so it unbalances the parentheses.

Your erroneous result seems to be related to a bad simplification of factorials. After evaluating at x=1 and r=1, your summand has the subexpression

ex:= (1+k)!/k!^2/(2-k)!:

which simplifies as

simplify(ex);

Obviously this is wrong for any integer k. A workaround is assume(k::nonnegint).

simplify(ex) assuming k::nonnegint;

To get rid of all of the RootOfs at once, you could do

evalindets(newpar, RootOf, x-> allvalues(x)[1]);

First method (my preference):  second:= eval(t, test[2]);

Second method: second:= rhs(test[2][]) or rhs(test[2][1]) or op([2,1,2], test)

Unfortunately, there is no analog to the first method for selecting the first value; it only works on expressions of the form (variable = value).

The trick is to first use command combine to bring the A into the sum, and then use algsubs(A*W[n](t)=5, ...). Finally, use expand to bring constants back out of the sums.

A*B*Sum(W[n](t)*q[n](z), n= 1..infinity);
                     /infinity               \
                     | -----                 |
                     |  \                    |
                     |   )                   |
                 A B |  /     W[n](t) q[n](z)|
                     | -----                 |
                     \ n = 1                 /
combine(%);
                  infinity                   
                   -----                     
                    \                        
                     )                       
                    /     A B W[n](t) q[n](z)
                   -----                     
                   n = 1                     
algsubs(A*W[n](t)=5, %);
                      infinity           
                       -----             
                        \                
                         )               
                        /     5 B q[n](z)
                       -----             
                       n = 1             
expand(%);
                         /infinity       \
                         | -----         |
                         |  \            |
                         |   )           |
                     5 B |  /     q[n](z)|
                         | -----         |
                         \ n = 1         /

Let me know if that's enough explanation, or if you need a more-detailed example.

The default solution method rkf45 fails miserably at this task. (I used x(4)=10, u(4)=10, z(4)=10.) Using method= dverk78 gives about 5 digits of accuracy. I am trying some other methods.

restart:
sys:= {diff(u(t),t) = 4*u(t), diff(x(t),t) = 6*x(t)+2, diff(z(t),t) = x(t)+u(t)}:
Sol1:= dsolve(sys union {u(4)=10, x(4)=10, z(4)=10}, numeric):
Sol1(0);
         [t = 0., u(t) = HFloat(1.11364764942213e-6), 
           x(t) = HFloat(-0.3333333330052356), 
           z(t) = HFloat(7.1111113895777045)]

Sol2:= dsolve(sys union eval({u(0)=u(t), x(0)=x(t), z(0)=z(t)}, Sol1(0)), numeric):
Sol2(4);
          [t = 4., u(t) = HFloat(9.742963826802734), 
            x(t) = HFloat(7.58924709847216), 
            z(t) = HFloat(9.533948806446041)]

Sol3:= dsolve(sys union {u(4)=10, x(4)=10, z(4)=10}, numeric, method= dverk78):
Sol3(0);
     [t = HFloat(0.0), u(t) = HFloat(1.12536549286537e-6), 
       x(t) = HFloat(-0.33333333294284945), 
       z(t) = HFloat(7.111111392517568)]

Sol4:= dsolve(

     sys union eval({u(0)=u(t), x(0)=x(t), z(0)=z(t)}, Sol3(0)),

     numeric, method= dverk78

):
Sol4(4);
      [t = HFloat(4.0), u(t) = HFloat(9.999698295626974), 
        x(t) = HFloat(10.000134201687256), 
        z(t) = HFloat(9.999946940854626)]

In terms of elegance and simplicity of the algorithm used, it'd be hard to beat rem:

f:= x^5+a*x^4+b*x^3+c*x^2+d*x+e:
rem(f, x^3-1, x);
                          2                    
                 (c + 1) x  + (a + d) x + b + e

@casperyc 

It really surprises me that this would happen because I simplified it.

Here's a simple example that shows why not to simplify. These are the same type of expressions as occur in your Matrices (the ones in your Matrices are longer).

simplify(1+(a+b+c+d)*(e+f+g+h)*(i+j+k+l));

 a e i + a e j + a e k + a e l + a f i + a f j + a f k + a f l
 + a g i + a g j + a g k + a g l + a h i + a h j + a h k
 + a h l + b e i + b e j + b e k + b e l + b f i + b f j
 + b f k + b f l + b g i + b g j + b g k + b g l + b h i
 + b h j + b h k + b h l + c e i + c e j + c e k + c e l
 + c f i + c f j + c f k + c f l + c g i + c g j + c g k
 + c g l + c h i + c h j + c h k + c h l + d e i + d e j
 + d e k + d e l + d f i + d f j + d f k + d f l + d g i
 + d g j + d g k + d g l + d h i + d h j + d h k + d h l + 1

Before simplification, there are 2 multiplications and 10 additions; after, there are 128 multiplications and 64 additions. So the complexity increased by a factor of 16. One way to guard against this is to include the size modifier to simplify.

simplify(1+(a+b+c+d)*(e+f+g+h)*(i+j+k+l), size);
      1 + (a + b + c + d) (e + f + g + h) (i + j + k + l)

But, based on an extremely brief viewing of some of the entries in your Matrices, I guess that simplify will provide no benefit in your case, even with the size option.

The functional way:

f:= x-> exp(-exp(x^2)):
plot(z-> Int(f, 0..z), -3..3);

Note the capital I in Int and the lack of a variable appearing before the ranges.

I assume that you intended for g to be considered a function of x, y, z. In that case, you should change Diff(g, x) to D(g)(x) or Diff(g(x,y,z), x), and likewise for y and z.

When you simplify the expressions in these matrices, they become incredibly long---essentially, they are expanded as polynomials in the ~ 10-20 variables. The only change I made was to remove the line kappa:= simplify(kappa). With this change, I was able to eval the K=8 case (the 254x20) in 0.25 seconds, and I did the K=10 (1022 x 24) case in 2.5 seconds and the K=11 (2046 x 26) in 8.5 seconds. The mykappa itself also runs much, much faster without the simplify. So, there, I gave you a factor of about 1000 improvement.

Not simplifying should also give you greater accuracy in the floating-point evaluation---fewer operations means less round-off error.

With some research in Wikipedia, I found an algorithm for counting the Eulerian circuits (not paths) in directed graphs (digraphs). It's called the BEST algorithm, BEST being the initials of the four authors' surnames. A digraph has an Eulerian circuit iff for every vertex v, indegree(v) = outdegree(v). When these are equal, we will just call it deg(v). Denote by N(G) the number of spanning trees of G. Then the number of Eulerian circuits is

N(G)*mul((deg(v)-1)!, v in G)

(For the graph in question, clearly all the factorials are 1.) Maple has a command for counting the spanning trees: GraphTheory:-NumberOfSpanningTrees. Applying the BEST algorithm in this case, we get 

restart:
macro(GT= GraphTheory):
G1:= GT:-Graph( {[C,G], [G,C], [T,G], [G,T], [A,T], [C,A], [T,C]}):
GT:-NumberOfSpanningTrees(G1)*mul((v/2-1)!, v in GT:-DegreeSequence(G1));
                               3

So, there are three Eulerian circuits.

(The division by 2 in the factorial is because DegreeSequence adds together the InDegree and OutDegree.)