HVNorm:= A-> VectorCalculus:-Norm(A[..-2])

The way that you constructed S4, it's vertex 1 in the center with degree 3, not vertex 4 like you show.

Do you understand, mathematically, what it means for a group to be generated by certain of its elements? The specified group is all 6 permutations of (2, 3, 4). It can be generated by the transpositions (2, 3) and (3, 4).

Object modules are extensively used by packages that have been written since objects were introduced. These include:

Calendar, ColorTools, DataFrame, Finance, GraphTheory, GroupTheory, Iterator, MutableSet, PolyhedralSets, PolynomialIdeals, TimeSeriesAnalysis.

It's the dominant paradigm now, and surely some programmmers at Maplesoft wish they could go back and rewrite some older packages in object form. IMO, Statistics sorely needs an OO overhaul. And any package that overloads infix operators could be improved by using objects.

It can be done with coeffs (not coeff) like this:

Coeffs:= proc(P, V)
local C,t;
    C:= coeffs(collect(P, V, 'distributed'), V, t);
    table(sparse, [t]=~ [C])
end proc
: 
C:= Coeffs(ODE, {Q(x), diff(Q(x),x)}):
for i from 0 to 1 do for j to 6 do c[i,j]:= C[diff(Q(x),x)^i*Q(x)^j] od od;
eval(c);

I don't know what you want as a final result. Here I give a minimal acceptable solution of the problem. My code produces a 10-second trajectory (of the bob only, not the rod or chain) for some initial conditions that I made up. Once you know the basic syntax of plotting commands and numeric solution of IVPs, the only tricky part to this is converting the solution function from spherical to cartesian coordinates. I put a comment over that part.

Animations are also possible.

params:= [ 
    l=1,    #length or radius (m)
    g= 9.8,    #gravity (m/s^2)
    theta0= Pi/3, phi0= 0,   #initial position (radians)
    thetap0=.1, phip0=1   #initial angular velocities (radians/s)
]:
eqs:= 
    diff(theta(t),t$2) = diff(phi(t),t)^2*sin(2*theta(t))/2 - g*sin(theta(t))/l,
    diff(phi(t),t$2) = -2*diff(theta(t),t)*diff(phi(t),t)*cot(theta(t))
;
ics:= phi(0)=phi0, D(phi)(0)=phip0, theta(0)=theta0, D(theta)(0)=thetap0:
sol:= dsolve(eval({eqs, ics}, params), numeric):
plots:-odeplot(
    sol, 
                                            #____the tricky part___#
    eval(changecoords([x,y,z]$2, spherical, [l, phi(t), Pi-theta(t)]), params), 
    t= 0..10, numpoints= 10^4, axes= normal, labels= [``$3],
    scaling= constrained
);

Pay particular attention to how derivatives are entered in differential equations---e.g., diff(phi(t), t$2)---and how they're entered in initial conditions---e.g., D(phi)(0) = ....

The implicitplot command is supposed to handle this (the system is polynomial), but I can't get it to work. I get an empty plot no matter what I tried with implicitplot. So, I resort to fsolve and regular plot. This gives 1 solution of x and y for each t although there may be multiple solutions. By adjusting the ranges, you can get others.

XY:= proc(T) 
option remember; 
    fsolve(eval({eq1,eq2}, t= T), {y= -1..0, x= -0.004..0})  
end proc:
X:= proc(T) try eval(x, XY(T)) catch: undefined end try end proc:
Y:= proc(T) try eval(y, XY(T)) catch: undefined end try end proc:
plot([X, Y, 0.1..2], axes= boxed, labels= [x,y]);
​​​​​​

To get other solutions, adjust the x range (-0.004..0), the y range (-1..0), and the t range (0.1..2).

The function is periodic with half period P:

P:= 33*Pi/1000/sqrt(1122):

The overall RMS of a periodic function equals the RMS over one period:

RMS:= evalf(sqrt(1/2/P*Int(x(t)^2, t= -P..P)));
                       0.000003200233260

See the Wikipedia article "Root mean square".

Here are problems 1, 3, and 4 plotted in parametric spherical coordinates:

x = rho*cos(theta)*sin(phi), y = rho*sin(theta)*sin(phi), z= rho*cos(phi).

So, rho is the distance from the point to the origin, theta is a longitudinal angle, and phi is the angle that the position vector of the point (i.e., the vector from the origin to the point) makes with the positive z-axis. In all cases, these are given as a list of three values [rho, theta, phi] (in that order!) where one of three parameters has been eliminated and replaced either with a constant or an expression depending on the remaining two parameters.

1. The region between the cone z = sqrt(x^2+y^2) and the plane z = 4:

plots:-display(
    #lateral surface of cone:
    plot3d(
        [rho, theta, Pi/4], rho= 0..4*sec(Pi/4), theta= -Pi..Pi,
        coords= spherical
    ),
    #circular end cap:
    plot3d(
        [4*sec(phi), theta, phi], theta= -Pi..Pi, phi= 0..Pi/4, 
        coords= spherical
    )
);

3. The region bounded by the hemisphere z = sqrt(25-x^2-y^2), the cylinder x^2 + y^2 = 4, and the plane z = 0.

plots:-display(
    #hemisphere:
    plot3d(
        [5, theta, phi], theta= -Pi..Pi, phi= arcsin(2/5)..Pi/2,
        coords= spherical
    ),
    #cylinder:
    plot3d(
        [2*csc(phi), theta, phi], theta= -Pi..Pi, phi= arcsin(2/5)..Pi/2,
        coords= spherical
    ),
    #annular end cap in xy-plane:
    plot3d(
        [rho, theta, Pi/2], rho= 2..5, theta= -Pi..Pi, 
        coords= spherical
    ),
    scaling= constrained
);

4. The region bounded by the hemisphere z = sqrt(16-x^2-y^2) and the paraboloid z = x^2 + y^2:

phi_int:= eval(phi, solve({cot(phi)*csc(phi)=4, phi>0, phi<Pi/2}));
plots:-display(
    #paraboloid:
    plot3d(
        [cot(phi)*csc(phi), theta, phi], theta= -Pi..Pi, phi= phi_int..Pi/2,
        coords= spherical
    ),
    #sphere:
    plot3d(
        [4, theta, phi], theta= -Pi..Pi, phi= 0..phi_int, coords= spherical
    ),
    scaling= constrained
);

I think that this functionality is covered by the command Units:-AddBaseUnit, but I don't have experience with this.

When applying evalf to the complicated algebraic expressions that usually come from the exact solutions of high-degree polynomials with parameters, one often gets small imaginary parts. If the expression is known to be real, these small parts can be considered to be 0. An example is

[solve(x^3 = a, x)];
eval(%, a= -1);
evalf(%);
   [                                                  -10                                 ]
   [0.5000000001 + 0.8660254037 I, -1.000000000 + 3 10    I, 0.5000000000 - 0.8660254039 I]

The imaginary part of the middle value should be 0.

The plot command is already taking this into consideration.

Now, if you still want to treat these values as truly imaginary, let me know.
 

For an unknown reason, this 10-year-old thread has been refreshed. So, I thought that this would be a good place to point out that the Question can now be Answered with a stock solution built-in to sort:

sort(x, 'output'= 'permutation');
             [1, 3, 2, 4]

Another way of thinking of it is that we want the 8th roots of 1+I:

(M,A):= (abs, argument)(1+I):
plot(((M^(1/8)*~[cos,sin])~([seq]((2*Pi*k+A)/8, k= 0..8))), style= pointline);

Maple does not automatically recognize vectors (or matrices) with equal entries as being equal, but they will recognized if they are converted to lists. If at the end you do

K1:= [seq]~(x, x=~ K);

then they'll be converted to a set of 2-lists. If you do

K1:= `<,>`~([seq]~(x, x=~ K)); 

then they'll be converted to a set of 2-vectors.

Are you trying to find the orbit of L under mulplication by M?

Your major problem is that doSum doesn't use its parameter first.

I suspect that this is an issue related to installation and licensing. I've seen numerous similar reports where that was the issue. You have entered the differential equation correctly. 

Try going to menu Tools -> Check for Updates.