dsolve(ode)  finds it. It is given implicitely.
You may use dsolve(ode, explicit)  but an inherent RootOf will appear. 

What you wrote is far from being the roots of the equation 2*x - tan(x) = 0.
This equation has infinitely many roots (just plot f to see this).

x=1/2*(tan(x))  is an equivalent form of the initial equation.

x=arctan(2x)  is also equivalent but only for -Pi/2 < x < Pi/2;  it has only three roots.

The use of unapply is the best choice.
The first construct does not work because the body of a procedure cannot be altered (without special tools).
A workaround, if you insist in using the -> construct:
seq( (a -> (x -> a*x))(a), a=1..3);

The option view=[...] works as documented. The fact that view=ymin..ymax   also works should be considered as a bonus (and not blamed). Actually for plot3d this is mentioned in the help.

Of course a more organized help system would be nice. In the past, searching the help system was easier and more relevant. It often happens to be more userful a Google search!

- Compute Rank(J) to see whether J is rank defficient.

- If you want to generate matrices with a given rank, multiply a diagonal matrix by invertible matrices. E.g.

with(LinearAlgebra):
A,B:='RandomMatrix(6,6)'$2;
R:=DiagonalMatrix([1$4],6);  #rank=4
J:=A.R.B;
Rank(%);


[A or B could be singular (improbable); then the rank could be <4].

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.