@STHence The worksheet is your worksheet with a few lines commented out and with my additions, which are in red Maple notation.

MaplePrimes13-03-19D.mw

The computation of the integral with the given relative tolerance (epsilon) takes 1-2 minutes (I didn't time it), so if you want a graph of G(t) (or sqrt(G(t))) you may want to use the plot options adaptive = false and numpoints = whatever you have the patience for. I didn't try that either.

Update: An idea for plotting G(t) or sqrt(G(t)) which avoids repetitive calculation of the same integral is to begin by deciding which t-values you want to use in plotting. Say 0, t1, t2, ..., tN.
Then find the integrals int( g ,t0..t1), int( g ,t1..t2),  ....
Below I have done that with an equidistant sequence of t's.

t0:=time():
N:=15: #Modest
for i from 1 to N do
   t1:=2*(i-1)/N;
   t2:=2*i/N;
   val[i]:=evalf(Int(g,t1..t2,epsilon=1e-3));
end do;
time()-t0; #10-11 min on my not very fast machine
plot([[0,0],seq([2*i/N,add(val[j],j=1..i)],i=1..N)]);
plot([[0,0],seq([2*i/N,sqrt(add(val[j],j=1..i))],i=1..N)]); #sqrt(G(t))

@STHence The worksheet is your worksheet with a few lines commented out and with my additions, which are in red Maple notation.

MaplePrimes13-03-19D.mw

The computation of the integral with the given relative tolerance (epsilon) takes 1-2 minutes (I didn't time it), so if you want a graph of G(t) (or sqrt(G(t))) you may want to use the plot options adaptive = false and numpoints = whatever you have the patience for. I didn't try that either.

Update: An idea for plotting G(t) or sqrt(G(t)) which avoids repetitive calculation of the same integral is to begin by deciding which t-values you want to use in plotting. Say 0, t1, t2, ..., tN.
Then find the integrals int( g ,t0..t1), int( g ,t1..t2),  ....
Below I have done that with an equidistant sequence of t's.

t0:=time():
N:=15: #Modest
for i from 1 to N do
   t1:=2*(i-1)/N;
   t2:=2*i/N;
   val[i]:=evalf(Int(g,t1..t2,epsilon=1e-3));
end do;
time()-t0; #10-11 min on my not very fast machine
plot([[0,0],seq([2*i/N,add(val[j],j=1..i)],i=1..N)]);
plot([[0,0],seq([2*i/N,sqrt(add(val[j],j=1..i))],i=1..N)]); #sqrt(G(t))

Thanks!
After I posted the answer I got to think (my usual order of doing things). Why 10 when 1*3 + 3*2 = 9?

Thanks!
After I posted the answer I got to think (my usual order of doing things). Why 10 when 1*3 + 3*2 = 9?

@Markiyan Hirnyk Yes I see that my integrand f was wrong.

Now the question is which was intended:

f:=1/sqrt(x^2+1)*ln(x^2+sqrt(x^2+1));

or

f:=1/sqrt(x^2+1)/ln(x^2+sqrt(x^2+1));

In the first case the same change of variables works. In the second it doesn't.
Unfortunately, I guess he meant the second.

@Markiyan Hirnyk Yes I see that my integrand f was wrong.

Now the question is which was intended:

f:=1/sqrt(x^2+1)*ln(x^2+sqrt(x^2+1));

or

f:=1/sqrt(x^2+1)/ln(x^2+sqrt(x^2+1));

In the first case the same change of variables works. In the second it doesn't.
Unfortunately, I guess he meant the second.

@adel-00 The problem of 3 real solutions for u as well as for v in a certain a-range must be dealt with.

Actually I think it is better to revert to using solve. So here is the whole thing from the start:

restart;
eq1:=(alpha+(l+alpha)*u+alpha*k*u^2)*a=
u*(alpha+(l+alpha)*u+alpha*k*u^2)*(1+l*alpha*b/((alpha+(l+alpha)*u+alpha*k*u^2)));
eq2:=a=(u*(1+l*v/(1+u+k*u^2)));
factor((rhs-lhs)(eq1));
eq1:=collect(%,u);
params:={l=1,alpha=1,b=100,k=50};
U:=[solve(eval(eq1,params),u)]: #3 solutions for u
plots:-complexplot(U,a=0..20,style=point); #plot in the complex u-plane
vua:=eval(solve(eq2,v),params); #v expressed in terms of u and a
V:=eval~(vua,u=~U): #the 3 solutions for v in terms of a
plots:-complexplot(V,a=0..20,style=point); #plot in the complex v-plane
#plot skips imaginary values of u, but notice 3 colors.
plot(U,a=0..20,labels=[a,u]);
## I'm assuming that this is what you were really after:
plot(V,a=0..20,labels=[a,v]);

There was a plotting problem in Maple 16.02, which may affect the plots above.
See the correction at
http://www.mapleprimes.com/questions/143131-DEPlot-Error-When-Plotting-Solution-Curve?submit=143189#comment143189

@adel-00 The problem of 3 real solutions for u as well as for v in a certain a-range must be dealt with.

Actually I think it is better to revert to using solve. So here is the whole thing from the start:

restart;
eq1:=(alpha+(l+alpha)*u+alpha*k*u^2)*a=
u*(alpha+(l+alpha)*u+alpha*k*u^2)*(1+l*alpha*b/((alpha+(l+alpha)*u+alpha*k*u^2)));
eq2:=a=(u*(1+l*v/(1+u+k*u^2)));
factor((rhs-lhs)(eq1));
eq1:=collect(%,u);
params:={l=1,alpha=1,b=100,k=50};
U:=[solve(eval(eq1,params),u)]: #3 solutions for u
plots:-complexplot(U,a=0..20,style=point); #plot in the complex u-plane
vua:=eval(solve(eq2,v),params); #v expressed in terms of u and a
V:=eval~(vua,u=~U): #the 3 solutions for v in terms of a
plots:-complexplot(V,a=0..20,style=point); #plot in the complex v-plane
#plot skips imaginary values of u, but notice 3 colors.
plot(U,a=0..20,labels=[a,u]);
## I'm assuming that this is what you were really after:
plot(V,a=0..20,labels=[a,v]);

There was a plotting problem in Maple 16.02, which may affect the plots above.
See the correction at
http://www.mapleprimes.com/questions/143131-DEPlot-Error-When-Plotting-Solution-Curve?submit=143189#comment143189

@adel-00 eq1 is a polynomial equation of degree 3. Maple gives 3 solutions when it realizes that fact and when they are real. That can complicate matters.
Here I find all three (complex) solutions for u and after that for v.

restart;
eq1:=(alpha+(l+alpha)*u+alpha*k*u^2)*a=
u*(alpha+(l+alpha)*u+alpha*k*u^2)*(1+l*alpha*b/((alpha+(l+alpha)*u+alpha*k*u^2)));
eq2:=a=(u*(1+l*v/(1+u+k*u^2)));
factor((rhs-lhs)(eq1));
eq1:=collect(%,u); #Now obiously a polynomial
params:={l=1,alpha=1,b=100,k=50};
#Finding 3 solutions for each a. I have made the increment in a 0.1 to see the changes more clearly.
cupoints:=[seq(fsolve(eval(eq1,params),u,complex),a=0..20,.1)]:
plots:-complexplot(cupoints,style=point);
#Making an animation:
Su:=seq(plots:-complexplot(cupoints[1..n],style=point,symbolsize=15),n=1..nops(cupoints)):
plots:-display(Su,insequence=true);
#Turning to the v-values (slightly repetitive):
cpoints:=[seq([[a,a,a],[fsolve(eval(eq1,params),u,complex)]],a=0..20,.1)]:
Vp:=eval(solve(eq2,v),params);
Sau:=[seq( zip~( (x,y)->[a=x,u=y],op(cpoints[i])),i=1..nops(cpoints))]:
cvpointsTemp:=ListTools:-Flatten([seq(eval~(Vp,Sau[i]),i=1..nops(cpoints))]):
#Removing disturbing Float(undefined):
cvpoints:=remove(type,cvpointsTemp,undefined):
plots:-complexplot(cvpoints,style=point);
#Animating the v-points
Sv:=seq(plots:-complexplot(cvpoints[1..n],style=point,symbolsize=15),n=1..nops(cvpoints)):
plots:-display(Sv,insequence=true);

Added: A simple graphical illustration showing when there are 1, 2, or 3 real solutions:

plots:-animate(plot,[(eval(eq1,params)),u=-1..7,-20..20],a=2.7..6.5);

@adel-00 eq1 is a polynomial equation of degree 3. Maple gives 3 solutions when it realizes that fact and when they are real. That can complicate matters.
Here I find all three (complex) solutions for u and after that for v.

restart;
eq1:=(alpha+(l+alpha)*u+alpha*k*u^2)*a=
u*(alpha+(l+alpha)*u+alpha*k*u^2)*(1+l*alpha*b/((alpha+(l+alpha)*u+alpha*k*u^2)));
eq2:=a=(u*(1+l*v/(1+u+k*u^2)));
factor((rhs-lhs)(eq1));
eq1:=collect(%,u); #Now obiously a polynomial
params:={l=1,alpha=1,b=100,k=50};
#Finding 3 solutions for each a. I have made the increment in a 0.1 to see the changes more clearly.
cupoints:=[seq(fsolve(eval(eq1,params),u,complex),a=0..20,.1)]:
plots:-complexplot(cupoints,style=point);
#Making an animation:
Su:=seq(plots:-complexplot(cupoints[1..n],style=point,symbolsize=15),n=1..nops(cupoints)):
plots:-display(Su,insequence=true);
#Turning to the v-values (slightly repetitive):
cpoints:=[seq([[a,a,a],[fsolve(eval(eq1,params),u,complex)]],a=0..20,.1)]:
Vp:=eval(solve(eq2,v),params);
Sau:=[seq( zip~( (x,y)->[a=x,u=y],op(cpoints[i])),i=1..nops(cpoints))]:
cvpointsTemp:=ListTools:-Flatten([seq(eval~(Vp,Sau[i]),i=1..nops(cpoints))]):
#Removing disturbing Float(undefined):
cvpoints:=remove(type,cvpointsTemp,undefined):
plots:-complexplot(cvpoints,style=point);
#Animating the v-points
Sv:=seq(plots:-complexplot(cvpoints[1..n],style=point,symbolsize=15),n=1..nops(cvpoints)):
plots:-display(Sv,insequence=true);

Added: A simple graphical illustration showing when there are 1, 2, or 3 real solutions:

plots:-animate(plot,[(eval(eq1,params)),u=-1..7,-20..20],a=2.7..6.5);

@adel-00 a=0 is a little special, so I left that out:

vpoints:=seq([a,subs(fsolve(eval({eq1,eq2},params),{u,v}),v)],a=1..20);
plot([vpoints]);
plot([vpoints],style=point);

# a = 0 is special because then the equations are satisfied if u = 0 and v is any number:

eval({eq1,eq2},params union {a=0,u=0});

#A slightly different approach:

V:=aa->subs(fsolve(eval({eq1,eq2},params union {a=aa}),{u,v}),v);
plot(V,0..20);

@adel-00 a=0 is a little special, so I left that out:

vpoints:=seq([a,subs(fsolve(eval({eq1,eq2},params),{u,v}),v)],a=1..20);
plot([vpoints]);
plot([vpoints],style=point);

# a = 0 is special because then the equations are satisfied if u = 0 and v is any number:

eval({eq1,eq2},params union {a=0,u=0});

#A slightly different approach:

V:=aa->subs(fsolve(eval({eq1,eq2},params union {a=aa}),{u,v}),v);
plot(V,0..20);

I tried this before I saw yours, but my version fails:

f:=proc(x,y) if not type([x,y],list(numeric)) then return 'procname(_passed)' end if;
      `if`(y>0,[x,4*y],[x,y])
end proc;

Edited: In the initial version of this comment I forgot 'return' . But the edited is the one intended and it doesn't work.

Edited once again: The version works if  _passed is replaced by x,y. Probably because _passed will be arguments passed to another procedure.

I tried this before I saw yours, but my version fails:

f:=proc(x,y) if not type([x,y],list(numeric)) then return 'procname(_passed)' end if;
      `if`(y>0,[x,4*y],[x,y])
end proc;

Edited: In the initial version of this comment I forgot 'return' . But the edited is the one intended and it doesn't work.

Edited once again: The version works if  _passed is replaced by x,y. Probably because _passed will be arguments passed to another procedure.

You say that you would like to find the analytic expression for that sum.

It may be that none such exists!