If beta is an integer, the system is polynomial.
But just look at the solution e.g. for beta=4:

solve(subs(beta=4,[f1,f2]),[x,y]);

With Maple there are better and easier to code methods.

restart;
Digits:=15:
ide:=x^2*(diff(y(x), x, x))+50*x*(diff(y(x), x))-35*y(x) - (1-exp(x+1))/(x+1)-(x^2+50*x-35)*exp(x)-int(exp(x*t)*y(t),t=0..1):
n:=10: nx:=16:
p:=x->add( a[k]*x^k,k=0..n):        #  edited
F:=unapply( eval(ide, [y=p]),x):     # (1)  edited
d:=evalf(add(F(k/nx)^2,k=1..nx)):
s:=Optimization:-NLPSolve(d,iterationlimit=20000):
yy:=eval(p(x),s[2]);   # approx solution

 numapprox[infnorm](eval(F(x),s[2]),x=0.1..1); #absolute error for ide
     0.676914366650072e-3

plot(yy,x=0..1);

(1)  Preben noticed an error here; it was essentially   F:=unapply( eval(ide, [y(x)=p(x),y(t)=p(x)]),x):

He also noticed that an exact solution of the ide is exp(x). Let's compare:

numapprox[infnorm](yy-exp(x),x=0..1);
       0.302777540774013e-3

I would write just two lines:

f:=D(F):
Int('f(x)',x=a..b)=int(f(x),x=a..b, continuous);

                    /b                       
                   |                         
                   |   f(x) dx = -F(a) + F(b)
                   |                         
                  /a                         

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);