The combine function works just fine for me:

> combine(2*sin((A + B)/2)*cos((A - B)/2));
                        sin(A) + sin(B)

The commutator in the tensor package takes vector fields (considered as rank one tensor fields) as input and produces a new vector field: it is the socalled Lie bracket of vector fields. It is therefore not applicable for matrix calculus as you envisage it. That being said, it is quite understandable that you may think that commutator() may be applicable for matrix calculus because, as far as I know, nowhere in Maple is something as common as the commutator of matrices to be found: that fact, if indeed true, never ceases to amaze me. So in many of my worksheets I just define my own commutator like

com := (A,B) -> A . B - B . A:

Then, for instance, you may test whether the Pauli matrices constitute a representation of the su(2) Lie algebra as follows:

with(LinearAlgebra):
pauli[1] := Matrix([[0,1],[1,0]]):
pauli[2] := Matrix([[0,-I],[I,0]]):
pauli[3] := Matrix([[1,0],[0,-1]]):
Equal(com(pauli[1],pauli[2]),2*I*pauli[3]),
Equal(com(pauli[2],pauli[3]),2*I*pauli[1]),
Equal(com(pauli[3],pauli[1]),2*I*pauli[2]);

You may even define your own jacobi function and test whether the Pauli matrices do obey the Jacobi identity as follows:

jacobi := (A,B,C) -> com(A,com(B,C)) + com(B,com(C,A)) + com(C,com(A,B)):
jacobi(pauli[1],pauli[2],pauli[3]);

You can do something like

eq := isolate(
	+(1/R1+1/R2)*(diff(V2(t), t))
	-(diff(V1(t), t))/R1+(V2(t)-V1(t))/L
	+C*(diff(V2(t), `$`(t, 2))) = 0
,V1(t)):
rhs(eq) = lhs(eq);

where I have broken the differential equation over several lines with the sole purpose to fit the screen in this forum.

But your use is incorrect. Try something like:

f := x -> x[1];
f(Vector([a,b,c,d]));

Even though being a theoretical physicist, I am also quite confused about the Physics package. For instance, take something as simple as the following code:

with(Physics):
Define(Faraday[mu,nu],antisymmetric={{mu,nu}}):
Faraday[1,1] := 10;

This I would expect to raise an error because Faraday, being antisymmetric, cannot be assigned any nonzero value at its diagonal elements. But there is no error raised. This is just one of several examples I do not understand. Probably, or rather, hopefully, it is only me using the package the wrong way. But only so, I suppose, because I have great difficulties getting to the bottom of the package by reading the help text.

Note added: Clearly my code below is not as compact as Douglas' beautifully compact code above is. I wrote it, though, because initially I could not figure out the use of his code (because from the help text I could not readily figure out the use of foldl): first after posting, naturally, I realized (good morning to me!) that for the L in his code I should substitute one of the lists of lists to be summed. Coming to realize it, my code is essentially the same as Douglas', only not as compact. I have come up with the following recursive code:

s := (L) -> `if`(nops(L) < 2,L[],
	s([[seq(seq(x + y,x = op(1,L)),y = op(2,L))],op(3..-1,L)])
)[]:

For your lists it correctly produces

L1 := [[2],[3]]:
L2 := [[4],[4],[5,1]]:
L3 := [[3,4,1],[2],[2],[6,7]]:
s(L1);
s(L2);
s(L3);
                               5
                             13, 9
                     13, 14, 11, 14, 15, 12

The values of x at the local minimum and local maximum, respectively, can be read of from the output of the following code:

f := (x) -> 9*x+21/(x-13):
x_extrema := [solve(diff(f(x),x) = 0,x)];
y_extrema := map(f,x_extrema);

Inspired by Robert Israels use of patmatch, which I did not know the existence of prior to his post, I have written the following code:

pairCombine := proc(x)
   local a,b,s,t;
   if patmatch(x,
      a::constant*cos(t::algebraic) + b::constant*sin(t::algebraic),'s'
   ) then
      map(assign,s);
      return sqrt(a^2 + b^2)*combine(
         'cos'(arctan(b,a))*cos(t) + 'sin'(arctan(b,a))*sin(t),trig
      )
   end if;
   x
end proc:
fullCombine := proc(x)
   local a,s,t,i,j,y;
   if type(x,`+`) then
      for i from 1 to nops([op(x)]) do
      for j from i+1 to nops([op(x)]) do
         y := pairCombine([op(x)][i] + [op(x)][j]);
         if
            patmatch(y,a::constant*cos(t::algebraic),'s') or
            patmatch(y,a::constant*sin(t::algebraic),'s')
         then return fullCombine(x + y - op(x)[i] - op(x)[j]) end if;
      end do;
      end do
   end if;
   x
end proc:

An example:

expr := +1/5*cos(s)-1/5*3^(1/2)*sin(s)-3/sqrt(2)*cos(t)-3/sqrt(2)*sin(t):
fullCombine(expr);

2/5*cos(1/3*Pi+s)-3*sin(1/4*Pi+t)

For the specific expression you give above you could do something like:

expr := 1/2*cos(t)-1/2*3^(1/2)*sin(t):
subs(t = t + Pi/3,expand(subs(t = t - Pi/3,expr)));

A different example where the shift in angle is Pi/4:

expr := 1/sqrt(2)*cos(t)-1/sqrt(2)*sin(t);
subs(t = t + Pi/4,expand(subs(t = t - Pi/4,expr)));

This probably is cheating, though, as it presupposes that you know the appropiate shift in angle.

Please disregard this post, the content of which I have deleted because it (on second thoughts) was nonsense; I misunderstood what you were asking for.

Most probably this is not answering your question, but I wonder about the formulation of the exercise itself: To me the piecewise position vector r is twice differentiable:

seq(LinearAlgebra:-Equal(
	map(limit,map(diff,p1,t$i),t=0,left),
	map(limit,map(diff,p2,t$i),t=0,right)
),i=1..3);

But maybe by "differentiable" is meant infinitely differentiable!?

Concerning output from dsolve: Writing

sol:= dsolve(des,ics);

in place of

dsolve(des union ics);

seems to solve the problem. PS: This comment was intended as a reply to Robert Israel, but became misplaced. I guess my coffee machine finished brewing a little too late :-).

To me there seems to be two problems with your plot p3: 1. Without really understanding what you want to do it seems that you need two nested seq's, not one; you cannot write something as

seq(something,(k,n)=(0..314,0..314));

which does raise the error "Error, unable to execute seq". Instead you have to write something like

seq(seq(something,k=0..314),n=0..314);

2. When isolating the pointplot part in p3 as (here broken differently to fit the screen)

pointplot({[
	3*sin(u)*cos(v)(k/50)(n/50),
	3*sin(u)*sin(v)(k/50)(n/50),
	3*cos(u)(k/50)(n/50)
],[
	2.9*sin(u)*cos(v)(k/50)(n/50),
	2.9*sin(u)*sin(v)(k/50)(n/50),
	2.9*cos(u)(k/50)(n/50)
]});

Maple raises the error: "Error, (in pointplot) incorrect first argument".

There is no need to involve the displayprecision. Just use evalf as illustrated by the following example:

Matrix(3,3,(i,j) -> evalf(sin(i*j/9),6));

When you copy-paste the output from this example to another application the number of decimal figures come out correctly. Note that the displayprecision only controls with how many figures floating point numbers should be displayed, not how many figures they are actually calculated to.

The anticommutator is defined as {A,B} = AB+BA. So {A,B} may be zero even though AB is nonzero. This happens if AB = -BA, i.e., when A and B anticommutes, as it is possible for fermions. Also, even though AB is nonzero it may produce a zero ket when acting on some ket, because AB is an operator which (analogous to a function) maps the domain space (all kets) to a disjunct direct sum of the null space (zero kets) and the range space (nonzero kets). The following relation holds: dim(domain space of AB) = dim(null space of AB) + dim(range space of AB). PS: Please note that I have no experience with the Physics package, so the above is written on quite general grounds.