@jud 

elem:=[Elements(g2)[]]:
c1:=CayleyTable(g1,elements=elem):
c2:=CayleyTable(g2,elements=elem):
add(abs(c1-c2));  #0

@jud It is possible, but is has not much sense. You may compare the Cayley tables. But you must be sure that the elements and their orders are the same, so, to be safe, the tables should be constracted by you.
C6, as a PermutationGroup, is generated by a 6-cycle in Symm(6). How could it be exactly equal to a subgroup in Symm(5)?

@ecterrab Thank you. However I consider the expression generated by D as very fragile and it needed your workaround.
The symbolic diff is mainly used for nonnegative integers j, k, so, why use GAMMA(-2*j) inside it?
Actually, the workaround does not work (division by 0) when j,k are assumed nonnegint.

restart
f := (x, y) -> 1/(2 + x*y^2);
ex1 := D[1 $ j, 2 $ k](f)(x, y)  assuming j::nonnegint, k::nonnegint;
ex2:=simplify(eval(ex1, [k = 4, Sum = add]), GAMMA);
eval(ex2, [j=2,x=0]);

@jud As I said, you must check for isomorphism, not "=".

restart;
with(GroupTheory):
AreIsomorphic~(SubgroupLattice(Symm(5), output=list), CyclicGroup(6));

@ecterrab I don't understand how D[1$j, 2$k](f)(x,y); is correct

The denominator is infinity for j::nonnegint  and  produces an Error for our j=2. The GAMMA at numerator too.

@acer Unfortunately we have a bug in D (for symbolic orders) here:

restart;
f := (x,y) -> 1/(2+x*y^2):
D[1$j, 2$k](f)(0,0):
eval(%,[j=2, k=4]):  value(%);    # 0  ???
D[1$2,2$4](f)(0,0);               # 6

@jud   ~  is used to construct elementwise operators (see ?elementwise) which are applied for each element of a container (i.e. list, set, rtable).
f~([a,b,c], u)   ==>   [f(a,u), f(b,u), f(c,u)]

@jud You are not right. G itself is a permutation group. You can see this by typing
lprint(G)

It has 5 generators. The same group can be defined with 3 generators:
GG := PermutationGroup(gen);
lprint(GG);

P.S. If you want to get rid of Perm, use convert:

convert(Perm([[1,2],[4,5,6]]), disjcyc);
        [[1, 2], [4, 5, 6]]

So,

convert~(gen, disjcyc);

@jud 
 

with(GroupTheory):
S5:=Symm(5):
L:=SubgroupLattice(S5, output=list):
map([IdentifySmallGroup],{L[]});
{[1, 1], [2, 1], [3, 1], [4, 1], [4, 2], [5, 1], [6, 1], [6, 2], [8, 3], [10, 1], [12, 3], [12, 4], [20, 3], [24, 12], [60, 5],   [120, 34]}

# [6,2] = C6, see: https://en.wikipedia.org/wiki/List_of_small_groups
 

@jud This simply means that this Mathematica command identifies (isomorphically) the subgroups.
The SubgroupLattice in Maple is also able to detect the subgroups. See also DrawSubgroupLattice.

@Carl Love  sys was copy-pasted from tomleslie's,  so he has typos too.

@Carl Love Nice, vote up.
dsolve can be made to work in a few seconds, but the last 2 solutions are expressed via Mathieuxxx functions.

restart;
sys:= [ 
[diff(v__0(t),t$2)+v__0(t)=0,  v__0(0)=1, D(v__0)(0)=0],
[diff(v__1(t),t$2)+v__1(t)=-v__0(t)^3,  v__1(0)=0,  D(v__1)(0)=0],
[diff(v__2(t),t$2)+v__2(t)=-v__0(t)^2*v__1(t),  v__2(0)=0,  D(v__2)(0)=0],
[diff(v__3(t),t$2)+v__3(t)=-3*v__0(t)^2*v__1(t)^2-3*v__0(t)^2*v__2(t),  v__3(0)=0, D(v__3)(0)=0],
[diff(v__4(t),t$2)+v__4(t)=-6*v__0(t)*v__1(t)*v__2(t)-3*v__0(t)^2*v__4(t)-v__1(t)^3,  v__4(0)=0,  D(v__4)(0)=0],
[diff(v__5(t),t$2)+v__5(t)=-6*v__0(t)*v__1(t)*v__3(t)-3*v__1(t)^2*v__2(t)-3*v__0(t)^2*v__4(t)-3*v__0(t)*v__2(t)^2,  v__5(0)=0, D(v__5)(0)=0]
]:
print~(sys);
s:=NULL:
for i to 6 do  s:=s, dsolve((eval(sys[i], [s]))); print([s][-1]) od:

@jalal The principal planes correspond to eigenvectors.
They are T[..,j] , j=1..3, (normalized).
But here there are infinitely many because there is a double eigenvalue.
I have included only the principal plane corresponding to the single eigenvalue. You may add other if you want.

restart;
with(LinearAlgebra): with(plots):
n:=3:
f:=-x^2 + 2*y^2 + 2*z^2 - 6*x + 4*x*y - 4*x*z - 8*y*z + 4*z - 12:
f:=eval(f, [x=x[1],y=x[2],z=x[3]]);
L:=10:
quad:=implicitplot3d(f, x[1]=-L..L, x[2]=-L..L, x[3]=-L..L, style=surface, scaling=constrained):
A:=VectorCalculus:-Hessian(1/2*f,[seq(x[i],i=1..n)]):
b:=eval(Vector( [ seq(diff(f,x[i]),i=1..n)]), [seq(x[i]=0,i=1..n)]):
c:=eval(f,[seq(x[i]=0,i=1..n)]):
X:=Vector([seq(x[i],i=1..n)]):
solve([ seq(diff(f,x[i]),i=1..n)],{seq(x[i],i=1..n)}); # the center
X0:= Vector[column]( eval([seq(x[i],i=1..n)],%) ):
J,Q:=Eigenvectors(A);
PJ:=sort(J, output=permutation);
T:=Matrix(GramSchmidt([seq( Q[..,PJ[j]],j=1..n)],normalized)): # T is orthogonal
fnew:=simplify( (T.X+X0)^+. A. (T.X+X0) + b.(T.X+X0) + c ); # ==> Hyperboloid of Two Sheets
col:=[red,yellow,blue]:
ax:=seq(arrow(X0, T[..,j], length=10, width=0.3, color=col[j]), j=1..n): 
p1:=implicitplot3d((X-X0)^+ . Q[..,PJ[3]], x[1]=-L..L, x[2]=-L..L, x[3]=-L..L, style=surface, scaling=constrained, transparency=0.5):
display(quad, ax, p1, orientation=[175,63,21], caption="Hyperboloid of Two Sheets");

P.S. If you are using the document mode (which I don't recommend for programming), click the ">_" button to create a prompt, then press F5 and paste the code there.

@C_R 

restart;
f:=arctan(y,x) + arctan(-y,x):
# exhaustive check
simplify(f) assuming x>=0,y::real;
simplify(f) assuming x<0,y>0;
simplify(f) assuming x<0,y<0;
simplify(eval(f,y=0)) assuming x::real;

                               0
                               0
                               0
                    (1 - signum(0, x, 1)) Pi

So, for reals, (f <> 0)  <==>   (x<0 and y=0)
 

@Cata 

S:=[seq]([vq[1], [solve(V__out/vq[1] = 1/sqrt((-m*x^2 + m + 1)^2 + (vq[2]*(x - 1/x))^2), x, useassumptions)]], vq in `[]`~(v__line, Q__s)) assuming x::positive;

plot([ [seq]([u[1],min(u[2])], u=S),  [seq]([u[1],max(u[2])], u=S) ]);