Enclose the solution in a list, then use map with subs.  Rather than explicitly enclosing the output of solve in a list, call solve with the variables in a list (rather than a set), that returns a list of list of solutions.  Thus

k := n*(phi+ln(-1/(-exp(phi)+exp(phi)*p-p)))/phi:
sols := solve(n*p*(1-p)=(sigma)^(2), [p]);  # note the use of [p] rather then {p}.
map(subs, sols, k);

You can sum all the columns at once with ArrayTools:-AllAlongDimension:

(**) test := Matrix(3, symbol=T);
                                       [T[1, 1]    T[1, 2]    T[1, 3]]
                                       [                             ]
                               test := [T[2, 1]    T[2, 2]    T[2, 3]]
                                       [                             ]
                                       [T[3, 1]    T[3, 2]    T[3, 3]]

(**) ArrayTools:-AddAlongDimension(test,2);
                                    [T[1, 1] + T[1, 2] + T[1, 3]]
                                    [                           ]
                                    [T[2, 1] + T[2, 2] + T[2, 3]]
                                    [                           ]
                                    [T[3, 1] + T[3, 2] + T[3, 3]]

Don't use determinants.  Have you tried LinearAlgebra:-Eigenvalues?

That returns a procedure which you can then evaluate for particular results:

dsys := { diff(x(t), t) = -4*y(t)*x(t)
          , diff(y(t), t) = -y(t)+2*y(t)^2-4*z(t)*y(t)
          , diff(z(t), t) = -z(t)-2*u(t)*z(t)
          , diff(u(t), t) = -u(t)+2*y(t)*u(t)-4*u(t)*z(t)+2*u(t)^2
        }:
ics := {x(0) = 1, y(0) = 1, z(0) = 1, u(0) = 1}:
integ := dsolve(dsys union ics, numeric):

integ(1);
   [t = 1., u(t) = 2.25032691557073, x(t) = 0.167727774585702,
             y(t) = 0.241445274528124, z(t) = 0.0394710068014504]

You can also use it to plot the result

plots:-odeplot(integ, [seq([t,f(t)], f in [x,y,z,u])], 0..1);

deqs := { diff(x(t),t) = 1/(1-x(t)), x(0) = 0 }:
integ := dsolve(deqs, numeric):
integ(1);
Error, (in integ) cannot evaluate the solution further right of .50000006,
probably a singularity
integ('last');
               [t = 0.500000060786008, x(t) = 1.00000008947060]

If sed is available do

sed -i 's/W//g' data 

Use add, rather than sum.  It provides an optional increment and doesn't require the forward quotes:

r := add(exp(-r*x)*l[x +1]*m[x+1], x = 0..50, 5);

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);