@gkokovidis It appears that in this case replacing the upper limit 5 by 5.0 doesn't actually give a result coming from numerical integration.
What happens is that numerical integration is indeed tried, but it fails after a valiant attempt, so it falls back on exact integration using method = ftoc (as it happens).
Try the following code and you will see this statement printed near the end:

" int: Numeric integration routines did not produce an answer, falling back to exact integration. "

restart;
infolevel[int]:=5:
g :=  sqrt(1/(1-x^2));
CodeTools:-Usage(int(g, x = 0 .. 5.0));

Using L:=int( g, x = 0..5)  followed by evalf(L)  is in this case much faster because time is not wasted trying to integrate numerically.
This fails because it doesn't fall back on exact integration (as designed).
 

restart;
g :=  sqrt(1/(1-x^2));
evalf(Int(g, x = 0 .. 5));
#int(f,x=0..5, numeric); # This is just an alternative syntax

I never tried those two buttons before. But they surely don't work for me either.
I'm using the exact same: Firefox 58.0.1 under Windows 10 (64 bit).

@acer Thanks for the explanation.
Just for the record, on my 64 bit Maple on Windows 10 I get the error
Error, (in Optimization:-LPSolve) invalid input parameter to external routine
for the 4 examples in your worksheet using hardware floats combined with method=activeset.

@tomleslie I wonder what your system is. I'm using Maple 2017.3 64 bit on Windows 10.
I checked the Maple versions 17.02, 18.02, 2015.2, 2016.2, and 2017.3 for problems with the examples in their respective help pages for LPSolve (which is being used by Maximize).
All but Maple 17.02 had critical failures (and at the end resulting in a crash) unless Digits:=16 (or above)!

Use
op(0,a);

But what is your purpose in doing
a:=0.312[1];

@tomleslie The usual meaning of linearity of the function f from one vector space to another is that f satisfies
(1) f(x+y) = f(x)+f(y) for all vectors x and y
(2) f(t*x) = t*f(x)  for all vectors x and scalars t.

Thus the expression 21+a*x considered as a function of a is not linear.

@AndrewG You probably meant:
" so it essentially does what Rouben Rostamian has done without using the weights option" .

@AndrewG Even using weights=<1,1,10^15,1,1,1> appears to do fine.

You need 8 boundary conditions since your system has two dependent variables of highest order 4.
Also you didn't give us the values of p and L.

@mmcdara You must have had a concrete positive integer n, e.g. 20.
It doesn't matter because what you see is just ToInert working on f(x,y,n) evaluated!
But when n=20 (e.g.) you get
f(x,y,n);
with output sqrt(x)*sqrt(y) so you just see the result of
ToInert(sqrt(x)*sqrt(y));
 

Your system is huge you say, but couldn't you give us a (very) simplified version containing the apparently problematic procedure proc(f(x,y), a, b, c, etc...)  (also simplified if possible).
In particular it would help if we knew the relation of the independent variable 't' to f(x,y) and a, b, c, ... ?
And what is the output of the procedure?

@Rouben Rostamian  Very nice.
Your arguments even hold when the derivatives of theta(r) and sigma(r) are just bounded in an interval (0,delta), for some delta > 0.

Not to my knowledge. In fact you basically answered the question yourself by saying:
" Since there are no general methods for constructing Lyapunov functions ... ".
But in doing the algebra involved in the construction of a Lyapunov function Maple is clearly a great help.
 

@Ramakrishnan The solution I gave will apply to the whole session (or all sessions in case of Apply Globally).
Personally I never used equation labels, so I use Apply Globally.

@Markiyan Hirnyk We should be aware that complexplot3d plots the absolute value:
 

plot(abs(sqrt(1/x)*erf(sqrt(x))), x = -2 .. 2, discont,color=red,thickness=3): p1:=%:
plots:-complexplot3d(sqrt(1/z)*erf(sqrt(z)), z = -2-2*I .. 2+2*I, grid = [60, 60]): p2:=%:
tr:=plottools:-transform((x,y)->[x,0,y]):
plots:-display(tr(p1),p2);