Using scanf is not the best choice.  A problem is that any characters that are not accepted remain in the queue, so the if a user typed several characters, then hit enter, typing either y or n on the next line won't work because there are previous characters that are being processed.  You might be able to get this to work using the extra format characters, say

scanf("%c%[^]")

however, I wasn't successful (that is, there were some cases that would cause problems).  Simpler is to use readline.  You might also limit the looping.  Say

askUser := proc()
local rep;
    to 3 do
        printf("Are you OK? (y/n): ");
        rep := readline('terminal');
        rep := rep[1];
        if rep="y" then return true; end if;
        if rep="n" then return false; end if;
    end do;
    return FAIL;
end proc:

You might consider the efficiency of the above technique.  Each time stack:-new is called, a new Maple table is allocated.  That takes time and space. Similarly, if you are using combinat:-randperm to randomize the initial order, then each call generates a new list.  The old data structures are eventually garbage collected, but that has a cost.  You might consider creating your own module for handling shuffling, one that reuses an existing structure.  Here is a simple mockup:

Cards := module()
export Shuffle
    ,  Deal
    ;

local deck, top, n;

    n := 52;

    deck := Array(1..n, i -> i);
    top := n;

    Deal := proc()
    local card;
        if top < 1 then
            error "deck is empty";
        end if;
        card := deck[top];
        top := top-1;
        return card;
    end proc;

    Shuffle := proc()
    local j,k;
        top := n;
        for j from n to 2 by -1 do
            k := RandomTools:-Generate(integer(range=1..j));
            (deck[j],deck[k]) := (deck[k],deck[j]);
        end do;
        NULL;
    end proc;

end module:


use Cards in
    Shuffle();
    printf("%a\n", [seq(Deal(), i=1..52)]);
end use;

A table is not a stack, though it can be used to emulate one.  Maple doesn't provide a built in stack data structure, but the procedures stack and SimpleStack allow you to create and emulate one.  For no clear reason, the SimpleStack constructor doesn't permit initializing the stack with a given contents.  You could just push the initial values onto the stack, or you can use the stack procedure, which permits initializing the stack.  For example:

(**) S := stack:-new(a,b,c):
(**) stack:-pop(S);
                                                        c

(**) stack:-pop(S);
                                                        b

(**) stack:-pop(S);
                                                        a

(**) stack:-pop(S);
Error, (in stack:-top) empty stack

If the Shuffle routine shuffles a list, then you can do

stack:-new(Shuffle(SList)[]):

Interesting.  Maple's expand procedure does not know about Im or Re. 

y := Im(a*t^2+b*t):
expand(y) assuming t::real;
                                      2
                                Im(a t  + b t)

That can be readily corrected by extending expand

`expand/Im` := proc(z)
    if z :: `+` then
        return map(Im, z);
    else
        return z;
    end if;
end proc:

Now

expand(y) assuming t::real;
                               2
                              t  Im(a) + t Im(b)

You can assign an equivalent procedure for handling Re.

Note, by the way, that you should not have to use eval in map(expand, eval(M)), unless M is of type matrix.  If so, then consider using the newer Matrix data structure, which supercedes matrix.

By color bar do you mean a visual mapping from colors to scalars?  That would only be possible if you used a shading scheme that was scalar (say zhue).  But I don't know that that is an available option, though it would be useful.  

As hirnyk points out, your piecewise's are probably wrong.  You probably want

piecewise(t<100, 0, t<=160, 1, 0)

Alternatively,

piecewise(t < 100 or 160 < t, 0, 1)

I don't know of a fast method.  For convenience, increment the parameters of interest in loop so that it is automated. Alas you probably cannot use the Threads package to do this in parallel because the procedures aren't certified thread safe.  If you have a multi-core machine consider firing up multiple Maple kernels, each handling a subset of the parameter values.

An interesting question is whether evalf/int is thread-safe.  I'll guess not, but if it were then you could first generate all the uu's for incremented parameter values, then integrate them in separate threads to compute S and B. 

One way around this, that avoids any checking, is to wrap the expression in a list:

for i to 5 do
   subs(x = a, [k||i]);
end do;

That avoids errors if k||i evaluates to NULL or an expression sequence.

You are generally better off using eval rather than subs, but the same issue arises.  That is, do

for i to 5 do
   eval([k||i], x=a);
end do;

You could do all at once with

eval([seq([k||i], i=1..5)], x=a);

It works here.  A better technique is to call solve with the variable enclosed in a set, that returns an expression sequence of sets of equations defining the solution (rather than an expression sequence of values).  That is, do

k1 := -(ln(q/(-1+q+exp(phi)))+phi)/phi:
eqn1 := (1-q)/q = mu:
sol := solve(eqn1,{q});
                                            1
                              sol := {q = ------}
                                          1 + mu

subs(sol,k1);
                                      1
                   ln(---------------------------------) + phi
                               /       1              \
                      (1 + mu) |-1 + ------ + exp(phi)|
                               \     1 + mu           /
                 - -------------------------------------------
                                       phi

You can use LinearAlgebra:-RowDimension (and ColumnDimension). Or LinearAlgebra:-Dimension for both.

Use simplify with side-relations.  Here is a simple procedure, Identify, to create the identifications:

Identify := proc( s :: set )
local i,j,n;
    n := nops(s);
    {seq(seq(s[i]=s[j], j=i+1..n), i=1..n-1)}
end proc:

p := (x+y+z)^3:
simplify(p, Identify({x^2*y, x^2*z, y^2*x}));

             2          2      2          2              3    3    3
          9 x  z + 3 x z  + 3 y  z + 3 y z  + 6 x y z + x  + y  + z

You could specify an additional variable to algsubs, that works if you are replacing n, however, you need to remove the sqrt,

algsubs((Variance(X3))=sigma^2,k3,n);

The Maple ^ operator is syntactically not associative (more precisely, non-associative), so x^y^z is ambiguous and raises an error.  In 2D math things are different in that the action of typing the caret (^) switches to superscript mode. The switch to superscript make it natural to make the operation right-associative.

Interpreting x^y^z as (x^y)^z is non-standard.  In 2D mode one can already get that effect by typing, in 2D, x^y*z, which creates x^(y*z), which is equivalent to (x^y)^z. 

I suspect that your interpretation comes from calculator usage, where typing x^y^z computes (x^y)^z.

Use printf:

printf("Stationary point %a is a %a\n", s[i], sptype);

You might consider using the inert Power function.  It is already known to Maple (used with mod).  Here is an extension to combine that does what you want

`combine/Power` := proc(expr1)
local pwrs, rest, p, L, b, e, expr;
    if expr1 :: `*` then
        expr := map(procname,expr1);
        (pwrs,rest) := selectremove(type,expr,'Power(anything,anything)');
        pwrs := `if`(pwrs :: `*`
                     , [op(pwrs)]
                     , [pwrs]
                    );
        pwrs := [ListTools:-Categorize((a,b) -> op(1,a)=op(1,b), pwrs)];
        return rest*mul('Power'(op([1,1],L), add(op(2,p), p in L)), L in pwrs);
    elif expr1 :: `+` then
        return map(procname, expr1);
    elif expr1 :: 'Power(anything,anything)^anything' then
        expr := map(procname, expr1);
        (b,e) := op(expr);
        return Power(op(1,b), op(2,b)*e);
    else
        return expr1;
    end if;
end proc:


combine(23*Power(a,b)^2 / Power(a,c) * Power(b,c));
                                       23*Power(a,2*b-c)*Power(b,c)