Try

roundcoeffs1(ggg, indets(ggg), 4);

The equation has an infinity of roots.
To find th n-th root (n>=0) you may use

nthroot := (n,R,L) ->  2*R/L*RootOf(tan(x)+tanh(x), x, -Pi/4+n*Pi):

nthroot(3,R,L): evalf(%);
     17.27875966*R/L

with(plots):
complexplot3d(z -> 1/(1-z), -1-I .. 1+I, axes=boxed);


This will plot |f(z)|  with color defined by argument(f(z))  in the square |x|,|y|<1.

If you want the unit disc |z|<1 only, use:

F := proc(z) local w; w := Re(z)*exp(Im(z)*I); 1/(1-w) end proc:
changecoords(complexplot3d(F, 0+0*I .. 0.95+2*Pi*I, axes=boxed), cylindrical);

Even without symmetries it is possible to find particular solutions.
pde:= diff(u(x,t),t)=alpha*diff(u(x,t),x$2)+f(t);
pdesol:=pdsolve(pde, u(x,t), build);

Then, for x=0, u(0,t)=f(t),  dsolve ==> f(t).

Finally, u(x,t) = a*x^2 + b*x - 2*alpha*a + exp(t)*d,  a,b,d being constants.

Your GL23 is represented as a permutation group.
It has two generators:  (1, 2)(3, 5)(4, 7)   and  (1, 3, 6)(2, 4, 8)
[written as product of disjoint cycles].

There is no matrix here.

with(ImageTools): with(plots):
m:=400: n:=600:    mm:=200: nn:=300:
p:=textplot([[1.6,1.6,"Why not?", 'font'=["times","roman",80]]], 'view'=[1..2, 1..2],axes=none,color=blue,size=[n,m]):
fn:=cat(kernelopts(homedir),"/whynot.png"):  fnr:=cat(kernelopts(homedir),"/whynot-rot.png"):
Export(fn,p):
A := Read(fn):
B:=Array(1..mm,1..nn,1..3, datatype=float[8],order=C_order,fill=1.0):

rot:=proc(a::float[8], A::Array(datatype=float[8],order=C_order),
                                     B::Array(datatype=float[8],order=C_order))
option autocompile;
local i::integer[4],j::integer[4],ii::integer[4],jj::integer[4];
local co:=evalhf(cos(a)*mm/m), si:=evalhf(sin(a)*nn/n);
for i to floor(m) do for j to floor(n) do
  ii:= round((i-m/2)*co-(j-n/2)*si +mm/2);
  jj:= round((i-m/2)*si+(j-n/2)*co +nn/2);
  if ii<=mm and ii>0 and jj<=nn and jj>0 then
     B[ii,jj,1]:=A[i,j,1];B[ii,jj,2]:=A[i,j,2];B[ii,jj,3]:=A[i,j,3] fi;
  od;od;
NULL;
end:
 
rot(Pi/6.,A,B);
Embed(B);
#Write(fnr, B);

I'd recommend first the same thing for the smaller group Symm(30). Here the generators are implemented and it has only 265252859812191058636308480000000 elements.

Try something like this:

A:=<a,b,0; 0,0,c>;

consts:=[a=2,b=3,c=4]:
A1:=eval(A, consts):
A2:=eval(A, c=66):
A1,A2;

simplify(expand(%));
collect(%,indets(%,name));  #optional

Among the methods used by fsolve there are:

For equations: Newton, Secant,  Dichotomic, inverse  parabolic interpolation.
For systems:  Newton, methods based on approximating the Jacobian and partial  substitutions.
Ostrowski's method (enhanced Newton) is also used sometimes.
If a method does not work, another one could be invoked.

In order to find what is effectively used, set
infolevel[fsolve]:=1;  #or greater if you want more details

Note that it is not difficult to find examples where fsolve fails but solve (or even the user, by hand) solves it.
E.g.
fsolve({(2*x+y+1)*exp(-(x+y-1)^2),   (3*x+2*y-1)*exp(-(x-y-1)^2)});

eq := z^2+y^3+x^4 - ln(x+y+z);
z__xy := implicitdiff(eq, z, x, y);
z__yx := implicitdiff(eq, z, y, x);
'z__yx' - 'z__xy' = z__yx - z__xy ;

with(plots):with(plottools):
display(disk([0,0], 1, color=red,transparency=0.5),disk([1,0], 1, color=green),scaling=constrained);

Probably you want:

with(plots):
ball := proc (a, b) fieldplot([x-a, y-b], x = -2 .. 2, y = -2 .. 2) end proc;
ball(1,1);

Your version makes no sense because in the procedure ball, x and y are is evaluated (substituted) with 1 producing
fieldplot([1, 1], 1 = -1 .. 1, 1 = -1 .. 1); #nonsense

We may use Poisson' formula:

uu:=x^2/2 + 1/(2*Pi)*(1 - (x^2+y^2)) *
      Int( 1/( (x-cos(t))^2+(y-sin(t))^2  )*(-cos(t)^2/2), t=0..2*Pi):
simplify(value(uu)) assuming x<1/2,y<1/2,x>-1/2,y>-1/2;

Of course, the result is valid for all (x,y) not only for |x|<1/2, |y|<1/2 but unfortunately the assume facility is not able to do this.
Note that value(uu) without assumptions ==> undefined ...        [???]

Note also that pdsolve is not able to find the correct arbitrary functions _Fnn

You should use lists instead of sets, otherwise the order of the elements is unpredictable.

f:=[1,2,3,4]:
h:=[1,2,4,5]:
select(i->f[i]=h[i], [seq(1..4)]);

maybe:

f[%];

or

select(i->i=h[i], f);