Do you mean something like this?

> with(Statistics):
   L := [3,1,4,1,5,9,2,6,5,3,5];
   mu:= Mean(L);
   sigma:= StandardDeviation(L);
   L1 := select(is, L, RealRange(mu-sigma, mu+sigma));
   X := RandomVariable(EmpiricalDistribution(Array(L1)));
   S := Sample(X, 1000);

The legend doesn't have to be a string: it can be a mathematical expression such as sin(omega*t).

I don't see any provision to change the font in a component label.  But via Unicode, all sorts of alphabets are included.  For example, using Character Map in Windows I was able to put Greek and Hebrew letters in a button label.  

The real problem, I think, that you want your procedure S to be passed the name X rather than the value of X.  Try this:

> S:= proc(Y::uneval) subs(K=Y,eval(T)) end proc;

If you would rather see alpha - beta instead of -beta + alpha,

> sort(alpha - beta, [alpha]);

alpha-beta

This sorts the expression in decreasing powers of alpha.  Now (in the remainder of this session) whenever alpha - beta appears by itself or as an operand of an expression, it will be given in that form rather than as -beta + alpha.

You might also try

> cos(beta - alpha) := 'cos(beta-alpha)';

after which cos(beta-alpha) will stay as cos(beta-alpha).

Something unexpected: in Maple 11.02 (Standard GUI, on Windows XP) Jason's code

> plots[textplot]([x, y, typeset(`#mover(mi( "V" ),mo("→"))`)]);

(with numerical values for x and y) produces the desired effect, but in Maple 12.01 (Standard GUI, Windows Vista) it leaves double quotes around the V.  Is this a
design change or a bug?

The question was why

 R := k-> (4*k^2-6*k+(-p-2)*(p-1))/(4*k^2-2*k) ;

has R(k) = 1 - 1/k + O(1/k^2) as k -> infinity.

It's really quite simple.  Divide both numerator and denominator by 
k^2:

R(k) = (4 - 6/k + (-p-2)*(p-1)/k^2)/(4 - 2/k)
     = (4 - 6/k + O(1/k^2))/(4 - 2/k)

Now 1/(4 - 2/k) = (1/4)/(1-1/(2*k)) = (1/4)*(1 + 1/(2*k) + O(1/k^2))
so
R(k) = (4 - 6/k + O(1/k^2))*(1/4)*(1+1/(2*k) + O(1/k^2))
     = 1 - 1/k + O(1/k^2)

Try this:

printVectorLabel := proc (x, y, VECTORNAME)
 plots[textplot]([x, y,
 typeset(`#mover(mi(` || VECTORNAME||  `),mo("→"))`)]);
end proc:

To plot a curve of the form f(x,y) = g(x,y), you can use implicitplot in the plots package.

I'll take Q = 10000.

> plots[implicitplot](18^2-y^2 = 10000*(ln((3000*(18-y))*sqrt(0.412e-4))-ln(sqrt(x)))/(0.412e-4*Pi),
   x=0..20, y = 17 .. 18,gridrefine=3,crossingrefine=3);

If you give pointplot3d a list of points and use the option style=line, it will connect those points in the order given.  In order to make the edges of  a cube, you could give it a list corresponding to a path traversing all edges of the cube.  Each point will be repeated at least twice.

I tried your differential equation with initial conditions c(0)=1, D(c)(0)=m for various values of m.  I found the result switched very rapidly from a solution going to -infinity to one that hits the singularity at c(x) = 4/3, near m = -84.2316172.  In fact, here are the solutions for m=-84.2316171 (in red) and m = -84.2316172 (in blue):

Due to roundoff error, I'm not sure how much to trust these results, but the switching point should be quite close to this.  For x near 0 (before the two solutions start to diverge noticeably), the solution to your problem should agree with these, and it should be very close to 0 for larger x.  Clearly the problem is extremely sensitive.

How many ways are there of choosing 3 teams out of 11?

How many ways are there of choosing 3 teams out of the 4 European teams?

You can try using simplify with a side condition.

> simplify((x+1)^6+2, {(x+1)^2 = a});

2+a^3

> simplify((x^2+2*x+1)^3+2,{(x+1)^2=a});

2+a^3

Assuming you have a linear PDE with one dependent variable and n independent variables, start with the n x n Matrix of coefficients A.  For example (if the PDE is in "D" form:

> pde:= x*D[1,1](u)(x,y) + y^2 * D[1,2](u)(x,y) + D[2,2](u)(x,y) + x*D[1](u)(x,y) + u(x,y) = x;

> A:= Matrix(2,2,(i,j) -> coeff((lhs-rhs)(pde),D[i,j](u)(x,y)));

Note that this "double-counts" the off-diagonal coefficients, so adjust it by doubling the diagonal entries.

> for i from 1 to 2 do A[i,i]:= 2*A[i,i] end do:

It may help to assume the variables are real.

> assume(x::real, y::real);

The PDE is elliptic if the matrix A is either positive definite or negative definite.  Of course (if the coefficients are not constant) this may depend on the variables.  IsDefinite should give you the condition needed for this to be true.

> with(LinearAlgebra):
> IsDefinite(A, query=positive_definite);

0 < x and 0 < 4*x-y^4

The 0 < x is redundant: but you can try

> solve(%, x);

{1/4*y^4 < x}

> IsDefinite(A, query=negative_definite);

x < 0 and -4*x+y^4 < 0

(which is never true)

> solve(%, x);

(returns nothing)

So we conclude that this PDE is elliptic for x > y^4/4.

The PDE is hyperbolic if A is indefinite.

> IsDefinite(A, query=indefinite);

not (0 <= 4*x-y^4 and 0 <= 2+2*x or 0 <= 4*x-y^4 and 0 <= -2-2*x)

> solve(%, x);

{x < 1/4*y^4}

The PDE is parabolic (I think: I've seen conflicting definitions) if A has rank n-1 and is (positive or negative) semidefinite.

Perhaps you were supposed to use 1/(1-cos(t)) and 3/(1+cos(t)) rather than what you wrote.  Those parentheses are important!  Otherwise the division takes precedence over the addition, so the result would be just 1-cos(t) and 3+cos(t).