This Answer is to provide a direct answer to your titular Question. I'm not posting this because I think that you should do permutations this way, but rather because passing and swapping pointers is particularly easy in Maple (especially Maple 2019 or later, as is coded below).

Procedure that swaps two pointers:

swap:= (a::evaln, b::evaln)-> ((a,b):= eval(b), eval(a)):
#Usage:
A:= 3:  B:= 7:
swap(A,B):
A,B;
                              7, 3

The evaln declaration (evaluate to a name) makes the procedure use the pointer rather than the value of the passed argument. The eval (evaluate) evaluates the pointer to its ultimate value. 

Procedure to swap (transpose) array elements:

ArraySwap:= (A::{Vector,Array}, i::integer, j::integer)-> 
    (A[[j,i]]:= A([[i,j]]))
:
#Usage:
A:= Array(1..9, k-> ithprime(k)):
ArraySwap(A, 3, 7):
[seq](A);
                [2, 3, 17, 7, 11, 13, 5, 19, 23]

Both of these procedures are simple enough that I'd likely just put the code inline instead of using explicit procedures. In particular, note that Maple's multiple-assignment operations do swaps without requiring any explicit intermediary variables.

You can include identical() in the typespec and initialize to () (both equivalent to NULL). Example:

P:= proc()
local a::{`=`, identical()}:= ();
    a
end proc:

Combinatorial objects can be generated one at a time with the Iterator package.

What purpose is served by using the inert form of Nabla in the composition? Why isn't (Nabla@@4)(f(x,y,z)) good enough for you? You can still make it inert by replacing f with %f.

Note that (Nabla@@4)(f(x,y,z)) = Divergence(Gradient(Divergence(Gradient(f(x,y,z))))) = 
((Divergence@Gradient)@@2)(f(x,y,z))

Hmm, this is a curiously simple Question for someone who's been on MaplePrimes for 6 years, so I hope that I'm not missing some subtlety. Anyway, 

P:= R*T/(V__m - b) - a/V__m/(V__m + b)/sqrt(T):
diff(P,T); diff(P, V__m $ 2);
 

Hmm. This is very bad, for the reasons that you metion. Apparently, the polymorphism doesn't extend to module names. There's a long-ago-deprecated module named process in your global namespace. You can do this:

my_car:-process(my_car);

I recommend that the object name be used as prefix in all cases simply to improve code readability if nothing else.

The fill value of the Array also needs to be of the specified type. The default fill value is 0, which is not of your specified type. So, do this

A:= Array(fill= sol, datatype= solution_class)

Download Arctangent_series.mw

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".