However, with that approach also

w();

and

w(77,bb,pp); #or whatever

would result in 99.

Could you give us text in MaplePrimes or upload a worksheet instead of a picture?

@Markiyan Hirnyk Nonlinearity has nothing to do with. The jacobian is the result of linearizing about the equilibrium point. There would be no point in linearizing if the system was linear.

I resent the tone of your comment, as I have done with some of your comments in the past to me and to other people.

@Markiyan Hirnyk Nonlinearity has nothing to do with. The jacobian is the result of linearizing about the equilibrium point. There would be no point in linearizing if the system was linear.

I resent the tone of your comment, as I have done with some of your comments in the past to me and to other people.

@Markiyan Hirnyk Stability questions are sometimes difficult to settle, yes. But very often it is a very simple matter as is the case with the two points found by fsolve (still assuming the order given in my answer):

subs(sol1,[v1,v2,u1]);
J(op(%));
LinearAlgebra:-Eigenvalues(%);
#So that point is asymptotically stable.
subs(sol2,[v1,v2,u1]);
J(op(%));
LinearAlgebra:-Eigenvalues(%);
#So that point is unstable.

@Markiyan Hirnyk Stability questions are sometimes difficult to settle, yes. But very often it is a very simple matter as is the case with the two points found by fsolve (still assuming the order given in my answer):

subs(sol1,[v1,v2,u1]);
J(op(%));
LinearAlgebra:-Eigenvalues(%);
#So that point is asymptotically stable.
subs(sol2,[v1,v2,u1]);
J(op(%));
LinearAlgebra:-Eigenvalues(%);
#So that point is unstable.

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