I am not quite sure what your objective is, but maybe the following procedure mySort (which uses the indexing procedure `index/makeIndex`) could be useful to you:

`index/makeIndex` := proc(Idx::list,Arr::Array,Ent::list)
	if nargs = 2 then Arr[op(Idx)]
	else Arr[op(Idx)] := (Idx = op(Ent))
	end if
end proc:
mySort := proc(M::Matrix)
   sort(
      ListTools:-Flatten(convert(Array(makeIndex,M),listlist)),
      (x,y) -> evalb(rhs(x) < rhs(y))
   )
end proc:

From a Matrix, mySort produces a sorted list of expressions in which the left hand sides are lists of the entrys row and column numbers, and the right hand sides are the entrys values. An example:

k := Matrix(3,3,(i,j) -> i*j + i^2 - j^2):
mySort(k);

If you assume that anything inside the trace that is not a Dirac gamma matrix is a commuting scalar (function), like p[mu] and k[nu] both are, then the following procedure should do the job:

myTrace := proc(x)
   if type(Physics:-Expand(x),`+`) then
      map(myTrace,Physics:-Expand(x))
   elif has(x,Dgamma) then
      remove(has,x,Dgamma)*Trace(select(has,x,Dgamma))
   else
      Trace(x)
   end if
end proc:

Try the following expressions (the first one being yours):

myTrace((p[mu]*Dgamma[mu]+m)*(k[nu]*Dgamma[nu]+m));
myTrace(
   (p[mu]*Dgamma[mu]+m)*
   (k[nu]*Dgamma[nu]+m)*
   (k[rho]*Dgamma[rho]+m)
);
myTrace(
   (p[mu]*Dgamma[mu]+m)*
   (k[nu]*Dgamma[nu]+m)*
   (k[rho]*Dgamma[rho]+m)*
   (k[sigma]*Dgamma[sigma]+m)
);

The issue here raised is the same as that raised in the thread Sequence of maps in which Robert Israel's post convinced me that using the function unapply is the appropiate solution:

a := array(1..3);
for i from 1 to 3 do
    a[i] := unapply(i*t,t)
end do;

Maybe the following is useful:

sol := {k = 1.,beta_star = 2.};
select(has,sol,k),select(has,sol,beta_star);
select(has,sol,k)[],select(has,sol,beta_star)[];
rhs(select(has,sol,k)[]),rhs(select(has,sol,beta_star)[]);

where lines two to four constitute three different levels of pulling out the appropiate parts of sol.

Note added: The suggestions below by Jay Treiman and Douglas Meade are better than this one of mine.

Maybe the following is useful:

isOdd := proc(f::algebraic,var::name)
	evalb(simplify(f + subs(var = -var,f)) = 0)
end proc:
isOdd(sin(x),x),isOdd(cos(x),x),
isOdd(sin(x),y),isOdd(cos(x),y);
                   true, false, false, false

For the solve part of your question, maybe the following help page would be useful:

?SolveTools,Inequality,LinearMultivariateSystem.

I am also unable to persuade Maple to properly evaluate the trace. I have tried to help Maple by explicitly telling it that p[mu] and k[mu] commute with the Dirac gamma matrices:

Setup(
   %Commutator(p[mu],Dgamma[nu]) = 0,
   %Commutator(k[mu],Dgamma[nu]) = 0
);

But to no avail, as Maple correctly cannot assume that they are scalars. The farthest I can get is using the following:

Trace(Physics:-Expand((p[mu]*Dgamma[mu]+m)*(k[nu]*Dgamma[nu]+m)));

But that, of course, is not satisfying. I have searched in vain for a way to tell Maple that p[mu] and k[mu] are commutative scalars, in the hope that then the Trace would let them pass outside it leaving behind only the Dirac gamma matrices. What is needed is to have Maple recognize that

Trace(p[mu]*k[nu]*Dgamma[mu]*Dgamma[nu]);

which emerges from the Expand-expression above when using the Setup above, belongs to the sixth possibility in the if-structure given at the help page ?Physics,Trace. But I do not know how to obtain that.

The help page ?trig clearly says that the arguments of the trigonometric functions must be given in radians.

Of course, if you just want to show your students the plots themselves you could do something like

f := 180 / Pi:
plot(sin(x / f),x = 0..360);
plot(cos(x / f),x = 0..360);

But I guess that is cheating, or what? If the students themselves should type the above, then, of course, they need to know about radians.

Maybe others here at mapleprimes have better ideas.

Concerning Exercise 4: You can define the following procedure:

Fouriercoeff := proc(f::algebraic,numTerms::posint)
   1/Pi*int(f,x=-Pi..Pi),
   seq([
      seq(1/Pi*int(f*F(k*x),x=-Pi..Pi),k=1..numTerms)
   ],F in [cos,sin])
end proc:

It is used as follows:

a0,a,b := Fouriercoeff(x^2,10);

Note that besides the function f, the procedure must be provided also a positive integer numTerms specifying the number of terms in the returned lists of Fourier coefficients.

Concerning Exercise 5: Using the above, define the expression:

expr := a0/2
	+ add(a[i]*cos(i*x),i=1..nops(a))
	+ add(b[i]*cos(i*x),i=1..nops(b));

From that expression it can be seen that the appropiate value is x = Pi, from which one can obtain for the famous Euler sum the value Pi^2/6.

I have written the following procedure which I hope that you will find helpful:

myFourier := proc(period,amps::list,numCoeffs::posint)
   module()
      export a0,aAreas,bAreas,aCoeffs,bCoeffs;
      a0 := `+`(amps[])/period;
      aAreas,bAreas := seq(Matrix(nops(amps),numCoeffs,
         (i,j) -> evalf(amps[i]*X(2*(i-1)*j*Pi/period))
      ),X in [cos,sin]);
      aCoeffs,bCoeffs := seq(Vector[row](numCoeffs,
         (i) -> `+`(convert(X[1..nops(amps),i],list)[])*2/period
      ),X in [aAreas,bAreas]);
   end module
end proc:

The numbers in your example are correctly reproduced writing

wave := myFourier(10,[14,18.7,9,4.1,6.7,6,6.3,8.4,4,2.9],5):
wave:-a0;
wave:-aAreas;
wave:-bAreas;
wave:-aCoeffs;
wave:-bCoeffs;

Some comments concerning formats: The number of coefficients is specified by numCoeffs. The amplitudes, as specified by amps, must be provided as a list. The areas wave:-aAreas; and wave:-bAreas; are returned as Matrices. The coefficients wave:-aCoeffs; and wave:-bCoeffs; are returned as row Vectors.

Does raising the number of digits being used for floating-point calculations, for instance Digits := 20; (the default being 10), remove the need for scaling up your unkowns?

To select the real roots you could write something like

select(is,[solve(2*x^5 - 5*x^4 + 3*x^3 - 3*x^2 + x + 2 = 0,x)],real);

To get rid of the RootOf you could write something like

map(allvalues,solve(
	{2*x^5 - 5*x^4 + 3*x^3 - 3*x^2 + x + 2 = 0,2*y - 3*x = 4},[x,y]
));

where the line is broken solely to fit the screen here at mapleprimes in 1024x768 resolution (which I am using).

According to the output of

interface(verboseproc = 2);
print(stats);

it seems that what you are asking for is

print(`stats/functions/describe`);

which leads to

print(`stats/find_function`);

which leads to

print(`stats/describe/functions/mean`);

for the mean.

You can copy-paste quite easily using either keyboard shortcuts or the mouse.