Except for D, none of the lettered problems have independent variables so plotting is pointless.  Maybe you are supposed to plot the partial sums?

A convenient way to locate the cause of the error is to use Maple's debugger.  Just wrap your code in a procedure, call stoperror("") to set the debugger to catch any error, then execute the procedure.  Here you could do

go := proc()
uses Statistics;
    N:=100:
    n:=100:
    X:=RandomVariable(Normal(500,50));
    for i to N do
        X||i:=Sample(X,n):
        mu||i:=1/n*sum(X||i(1,j),j=1..n):
        sigma||i:=sqrt(1/(n-1)*sum(((X||i(1,j)-mu||i))^(2),j=1..n)):
        Z||i:=(mu||i-500)*sqrt(n)/sigma||i:
    end do:
end proc:

stoperror(""):
go();

7     sigma || i := sqrt(1/(n-1)*sum((X || i(1,j)-mu || i)^2,j = 1 .. n));
(*DBG*) 

The debugger opens at the offending line.  To figure out what is happening, you can manually evaluate the expression by entering it into the debugger input. That is, entering

(*DBG*) sum((X || i(1,j)-mu || i)^2,j = 1 ..n))

generates the same error,

Error, unsupported type of index, j

so you can safely assume the problem lies there.  If aware of the sum/add distinction you could evaluate the expression using add and confirm the fix.  That is

(*DBG*) add(X||i(1,j),j = 1 .. n)
51468.2182490783

There is, alas, a bug in the procedure. I've submitted an SCR against it.

If f is truly a function, not just an expression, then you can use the Maple D operator to compute the derivatives.  While I probably wouldn't go this route (preferring to work with expressions), you can do

approx5f :=(f,a)-> f(a)+D(f)(a)*(x-a)+1/2*(D@@2)(f)(a)*(x-a)^2:
approx5f(sin, 3);
                                                             2
                 sin(3) + cos(3) (x - 3) - 1/2 sin(3) (x - 3)

evalf(%);
                                                                 2
            3.111097498 - 0.9899924966 x - 0.07056000405 (x - 3.)

While correct, this approach is a bit weird in that it accepts a function and returns an expression.  If you truly want to work with functions, it would make more sense to return a function (of x).  That could be done with

approx5f :=(f,a)-> unapply(f(a)+D(f)(a)*(x-a)+1/2*(D@@2)(f)(a)*(x-a)^2,x):
sin3 := approx5f(sin,3);
                                                                   2
          sin3 := x -> sin(3) + cos(3) (x - 3) - 1/2 sin(3) (x - 3)

sin3(1);                                                                
                                                   -sin(3) - 2 cos(3)

Is C a Vector or a list?  Regardless, you can use the same technique:

(**) C:=[x,y,a,z,b,c,d]:
(**) mask := [3,5,6,7]:
(**) C[mask];
                                 [a, b, c, d]

A simple approach is to split the polynomial, work on each part separately, then recombine:

(**) p := 1+(x+y)+(x+y)*x; 
                                                        p := 1 + x + y + (x + y) x

(**) (p1,p2) := selectremove(has, p, {x,y});
                                                      p1, p2 := x + y + (x + y) x, 1

(**) factor(p1)+p2;                         
                                                           (1 + x) (x + y) + 1

Move the space before the 'f' to after it. You are usually better off using seq rather than $ for indexed sequences.  Or use a range in an index, if applicable.  Or, use op. For example,

printf(cat("%12.5f "$4), seq(vij(1,1)[k], k=1..4));
printf(cat("%12.5f "$4), vij(1,1)[1..4][]);
printf(cat("%12.5f "$4), op(1..4, vij(1,1)));

If vij(1,1) evaluates to a Vector you could do

printf("%12.5f\n", vij(1,1));

Use a piecewise expression to represent f.  For example

f := piecewise( x < 1, x, x^2):
plot(f, x = 1..2);

Use the Maple exp function.  That is, exp(x) = e^x. I don't know what you mean by "showing e".

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