This little procedure will temporarily treat negatives as positives for any command, any expression. Specifically it does this:

1. For every negative n found, substitute -n*p, where the p is just a placeholder to indicate that a correction will be needed later;
2. Apply the command;
3. Change p to -1.

NegAsPos:= proc(C,e)
local `-1`;
   eval(C(subsindets(e, negative, `-1`*`-`), _rest), `-1`= -1)
end proc:

So, if e is your expression, do

NegAsPos(expand, e)
 

Ordinarily, the commands freeze/thaw or frontend are used for this sequence of operations:

1. Make a temporary substitution of placeholders for all subexpressions of a given type;
2. Do some command;
3. Replace the placeholders with what they were originally equal to.

That won't work in this case because the inner command (expand in this case) needs to know the absolute values of the negatives, which it wouldn't be able to see if they were hidden inside a placeholder.

The vertical-axis label can be moved up or down by post-padding or pre-padding (respectively) with newline characters `\n`, and it can be moved left by post-padding with spaces. I don't know a formula for the number of characters needed, but it would involve the font size (your 19), and the axis size (your 600). In this particular case, the maximum height of the y label is achieved by using 19 newlines. It is not possible to raise the label so that it's higher than the entire plot.

L:= 19:
plot(
   x, x= 0..2, labels= ['x', cat('y', `\n`$L)], 
   labelfont= [TIMES,19], size= [1200,600], axes= boxed
); 

As shown in the other Answers, it's possible to "fake" other plot components (such as title) so that it appears to be an axis label. In this manner, a fake label could be made to appear above the plot. All plot components that take strings (even tickmarks) can be padded with spaces or newlines. 

If the results differ in only the last decimal place (as for n=2 in your worksheet), then the issue is simply round-off error, as explained in NM's Answer.

On the other hand, if the results are very different (as for your cases n = 3, 4, 5), then the following explanation applies:

The expression cos(n*arccos(x)), for positive integer n, can be expanded into a polynomial of degree n, hence it has n roots. Each of solve and fsolve picks one of the roots in a different and somewhat arbitrary (although not random) manner. If you do 

p:= expand(cos(5*arccos(x)))

you'll get 16*x^5 - 20*x^3 + 5*x. If you apply fsolve to this, you'll get the 5 roots, because fsolve returns all real roots of (explicit) univariate polynomials, but only one root of other equations.

On the other hand, it's fairly easy to see that the exact solutions are cos(Pi/10), cos(3*Pi/10), cos(5*Pi/10), cos(7*Pi/10), and cos(9*Pi/10). I can get solve to give these answers, but it's a little bit tricky, and perhaps irrelevant to your coursework.

Your equation can be easily solved symbolically by hand in elementary functions using only the methods of 2nd-semester calculus. It has no relation to the equation "from the literature." Perhaps you are missing a factor?

ode:= r*diff(H(r),r$2)+diff(H(r),r)+(r*k^2-r^2*b^2/R^2) = 0:
dsolve({ode, H(R)=0, H(1/R)=R});
simplify(%) assuming k>0, R>0, b>0;

The solution is of the form P3(r) + C*ln(r), where P3 is a cubic polynomial and C is a constant (both of which depend on b, k, and R).

Eqs:= 
   ax^2+ay^2 = bx^2+`by`^2,
   bx^2+`by`^2 = cx^2+cy^2,
   (ax-bx)^2 + (ay-`by`)^2 = 3^2,
   (ax-cx)^2 + (ay-cy)^2 = 5^2,
   (bx-cx)^2 + (`by`-cy)^2 = 7^2
:
#Arbitrarily pick one to set to 0:
Sol:= solve(eval({Eqs}, cy=0));
SolA:= allvalues(Sol[1]);
plot(
   eval(
      [
         [[ax,ay], [bx,`by`], [cx,0], [ax,ay]],
         [abs(cx)*cos(t), abs(cx)*sin(t), t= -Pi..Pi]
      ],
      SolA[1]
   ),
   thickness= 4, scaling= constrained
);

You cannot control the order that multiplicands appear in a `*` term.

If instead you use the inert multiplication operator %*, then the order will remain as you specified. For example,

add(x[indx[]]  %* e[indx[]], indx= indxes);

If you view the output with 2-D output (by setting interface(prettyprint= 2)), then it'll look reasonably normal.

You asked:

• Am I correct in assuming that when the code I am using in the worksheet with is incorporated into a module in an mpl file and loaded, this will not be an issue?

That's absolutely NOT correct. This is a feature of the kernel, not of the GUI worksheet^[*1]. You cannot control the order, period.

[*1] Although there may be cases where the order displayed in the worksheet output differs from the order in which `*` objects are stored in the kernel. This I am not sure about. Nonetheless, you cannot control the order that they're stored in the kernel, which is what matters. You can always inspect this order with the op command. 

Grid:-Map, the second method you show, is the correct way. 

Try issuing this command first:

interface(prettyprint= 2);

However, this should've been set correctly by default, so I wonder how it got set incorrectly.

The lambda calculus is fully implemented in Maple through the arrow -> operator.  For example

(x-> y-> x^2+y^2)(3)(4)

Very good Question. Vote up.

The keyword _rest represents all the extra arguments that don't correspond to a parameter. If there are no extra arguments, then _rest = NULL.

f:= proc(x,y,a:=1) a*x*y end proc:
g:= q-> q(x, y, _rest):
g(f), g(f,3);
      x*y, 3*x*y

Note that it's necessary here to give a a default vaiue, so I chose 1.

You can write the expression by

diff(w(x,t), [x$2, t$alpha])

But I don't know how to solve a PDE containing such a derivative.

The error message indicates that option 'outputs' should be set to a list. So, change outputs= x to outputs= [x] or outputs= [x(t)]. 

I don't have MapleSim, so I can't test. 

It can be done by direct integration:

simplify(int(t^a*exp(-t^2-I*w*t), t= -infinity..infinity));

-(((((-1)^a*I - I)*cos(a*Pi/2) - sin(a*Pi/2)*((-1)^a + 1))*CylinderD(a, -sqrt(2)*w/2) - CylinderD(a, sqrt(2)*w/2)*(((-1)^a*I - I)*cos(a*Pi/2) + sin(a*Pi/2)*((-1)^a + 1)))*2^(-a/2)*exp(-w^2/8)*sqrt(Pi))/(4*sin(a*Pi/2)*cos(a*Pi/2))

If you want the transform of exp(-t^2)/t^a, just change a to -a. 

Is this what you want?:

restart:
co:= blue:
plots:-animate(
   plot, 
   [['seq'([[20-i,0],[0,i]], i= 0..k)], color= co, thickness= 3],
   k= 0..20, frames= 21, paraminfo= false
);

If instead you'd like the line segments to have constant length (for example, to simulate a ladder being pushed up a wall), that can be done.

The command

combinat:-randcomb(A,3) 

will draw a sample of size 3 without replacement from A. Kitonum's method draws a sample with replacement. Both methods are valid, but their probability distributions (hypergeometric and binomial, respectively) differ significantly when the population A is small. You'll need to decide which correctly models your situation.