A very useful feature of display is overrideoptions:
 

restart;
pt:=plottools:-line([1,2], [3,4],color=red,thickness=2,linestyle=dash); 
plots:-display(pt);
plots:-display(pt,overrideoptions,linestyle=1,thickness=5,color=blue);

restart;
with(IntegrationTools):
J1 := Int(h(x), x=-infinity..a) - Int(h(x), x=-infinity..b);
Ja:=Int(h(x), x=-infinity..a);
Jb:=Int(h(x), x=-infinity..b);
# Case 1: a>b.
Sab:=Split(Ja,b) assuming a>b;
Combine(Sab-Jb);
# Case 2: a<b. 
Sba:=Split(Jb,a) assuming a<b;
Combine(Ja-(Sba)) assuming a<b;
simplify(%);

MaplePrimes24-02-27_IntegrationTools_Combine.mw

NOTE: You don't need the assumption at the end, but it doesn't hurt.
You could add Flip(%); after simplify(%);
Then the results from both cases look the same.

restart;
###
with(LinearAlgebra):
###
v := <x, y>;
w:= <a,b>;
v^%T.v;
v^%T.w;
BilinearForm(v, v) assuming real;
##
## So you could e.g. do this:
M:=module() option package; export BilinearForm;
      BilinearForm:=proc(v::Vector,w::Vector) v^%T.w  end  proc;
end module;  
##             
with(M);
BilinearForm(<x, y>,<a,b>);
BilinearForm(<x, y,z>,<a,b,c>);
##
LinearAlgebra:-BilinearForm(<x, y>,<a,b>);

In the first procedure I would do this:
 

restart;
J1 := proc()
  local z:
  z := proc(u) if not u::realcons then return 'procname'(_passed) end if;
  fsolve(sqrt(x)=u, x) end proc:
  evalf(Int(z(u), u=0..1))
end proc:
##############
J1(); # 0.3333333333

In the next I would do likewise:

J2 := proc()
  local z:
   z := proc(q1, q2) if not [q1,q2]::list(realcons) then return 'procname'(_passed) end if;
     exp(
       2*(
         fsolve(1/2+(1/2)*erf((1/2)*x*sqrt(2)) = q1, x)
         *
         fsolve(1/2+(1/2)*erf((1/2)*x*sqrt(2)) = q2, x)
       )
       -
       fsolve(1/2+(1/2)*erf((1/2)*x*sqrt(2)) = sqrt(q1), x)^2
       -
       fsolve(1/2+(1/2)*erf((1/2)*x*sqrt(2)) = sqrt(q2), x)^2
    )
  end proc:
  evalf(Int(z(q1, q2), q1=0..1, q2=0..1))
end proc;

The execution of J2(), however takes an enormous amount of time (infinity?),
but fsolve has to work 4 times inside evalf/int for each value of q1=0..1, q2=0..1.

## Note:  If you comment out evalf(Int(z(q1, q2), q1=0..1, q2=0..1))  and replace it with
z(1/2,1/5) the execution works and returns 0.7300930564.
##
Note 2: I tried using narrower intervals:
evalf(Int(z(q1, q2), q1=0.1..0.9, q2=0.1..0.9))
J2() returned 0.4412675858 in a short time.

You are defining the procedure f to be the same as the procedure F so why not simply start with f:=F  ?

restart;
f := F;
((f@@(-1))@f)(x);
# Let us assume that y is defined this way

y := ((f@@(-1))@g)(x);
# When g is identical to f I would like to get y=x

'y' = eval(y, g=f); # y=x

A concrete example where F=ln:
 

restart;
f := ln;
((f@@(-1))@f)(x);
# Let us assume that y is defined this way

y := ((f@@(-1))@g)(x);
# When g is identical to f I would like to get y=x

'y' = eval(y, g=f); #y=x

I corrected a few errors, and here is what I have:
 

restart;
with(plots):
ode1 := diff(f(x), x$3) + 
        f(x)*diff(f(x), x, x) - lambda*diff(f(x), x) - diff(f(x), x)^2 + 2*beta*g(x) - Fr*diff(f(x), x)^2 
        = 0;
ode2 := diff(g(x), x$2) - diff(f(x), x)*g(x) + f(x)*diff(g(x), x) - 2*beta*diff(f(x), x) - lambda*g(x) - Fr*g(x)^2 
        = 0;
# You have 6 boundary conditions, but you can only have 5, unless one of the constants e.g lambda has to be determined.
# Your boundary conditions all have values 0, 
# which means that f(x) = 0, g(x)= 0 for all x on the interval is a (trivial) solution.
# To get a possible nontrivial solution you need to give an approximate solution:
# See ?dsolve,bvp
bcs := f(0) = 0, D(f)(0) = 0, f(5) = 0, D(f)(5) = 0, g(0) = 0, g(5) = 0;
a1 := [lambda = 0.2, beta = 0.2, Fr = 0.1];
####
sys1:=eval({ode1,ode2,bcs[1..5]},a1);# I arbitrarily removed the sixth bcs (g(5)=0).
res1:=dsolve(sys1,numeric);
odeplot(res1,[x,f(x)],thickness=5,view=-0.01..0.01);# As expected the zero solution.

To find a suitable approximate solution is not easy if you have no idea how f and g are expected to look.
You might have an idea?

Here is a proof that y(t) -> infinity as t -> infinity:
MaplePrimes24-02-07_odesys_infinity.mw

Here is the code directly:
 

restart;

sys_ode := diff(x(t), t) = -x(t)^2/(4*Pi*y(t)*(x(t)^2 + y(t)^2)), diff(y(t), t) = y(t)^2/(4*Pi*x(t)*(x(t)^2 + y(t)^2));
ics := x(0) = 1, y(0) = 1:

ode1,ode2:=sys_ode;

ode2/ode1;

ode:=diff(y(x),x)=-(y(x)/x)^3;

sol:=dsolve({ode,y(1)=1});

#From ode1 it follows that x(t) is decreasing from its initial value 1.
# sol shows that x cannot get lower than 1/sqrt(2), which is lower than 1. 
# Since x(t) > 1/sqrt(2), but decreasing we get by replacing x(t) by 1/sqrt(2):
ineq:=rhs(ode2) > subs(x(t)=1/sqrt(2),rhs(ode2));
ode3:=diff(y(t),t)=lhs(ineq);
sol3:=dsolve({ode3,y(0)=1});
limit(rhs(sol3),t=infinity); #infinity
# Thus y(t) from sys_ode will go to infinity too.

####### Graphical illustrations:
plot(rhs(sol),x=1/sqrt(2)..5,color=red,view=0..3); p1:=%: #Saving the plot

Digits:=15;
Sol := dsolve({sys_ode, ics}, {x(t),y(t)}, numeric,abserr=1e-16,relerr=1e-13,events=[[x(t)=5,halt]]);

Sol(100000);
1/sqrt(2.);
plots:-odeplot(Sol,[t,x(t)-1/sqrt(2)],99000..100000,caption=typeset("The difference between x(t) and ",1/sqrt(2)));
#Since x(t) is decreasing from 1 at t=0, you have to go backwards to get near to 1/sqrt(2)
plots:-odeplot(Sol,[x(t),y(t)], t=-43.1..100, color=blue,view=0..3,numpoints=10000);  p2:=%: # Saving the plot

plots:-display(p1,p2);# The last one is on top, so the color is blue
plots:-display(Array([p1,p2]));

Here is a simpler version. Let me know, if the result is what you expected:
 

restart;
with(LinearAlgebra):

xA := 1;
yA := 0;
xB := 0;
yB := 0;
xC := 0;
yC := 1;
Mat := Matrix(3, 3, [xA, xB, xC, yA, yB, yC, 1, 1, 1]);

Miv:=Mat^(-1); #Simpler
#Miv := MatrixInverse(Mat);

Miv^%T; #Simpler
#Transpose(Miv);
Miv^%T. <x, y, 1>;

phi := unapply(%,x,y);

for i to 6 do
    B || i := phi(xA || i, yA || i);
end do;

MaplePrimes24-02-05_LinearAlgebra.mw

Try this:

restart;

foo1:=overload([                   
                  proc(P1::list,P2::list,P3::list,$)
                      option overload(callseq_only);
                    [P1,P2,P3]
                  end proc,
                  proc(P1::list,P2::list,a::algebraic:=4,$)
                      option overload;
                    [P1,P2,a]
                  end proc
                       ]):
foo1([1,2],[3,4]); 
foo1([1,2],[3,4],x);
foo1([1,2],[3,4],[4,7]);

See ?overload where we read:
"If p1 raises an exception during the parameter type-checking phase only (not inside p1), and option overload(callseq_only) has been specified, then execution will proceed to p2(x,y)."

Note: See below for my comment on Ronan's original foo1.
Also try running my example above with option overload instead of option overload(callseq_only):
The results are the same. Thus the answer to the (revised) title is it didn't matter here.
Option overload(callseq_only) was introduced in Maple 16, but no example where it matters was given.

restart;
f :=piecewise(0 <= x and x <= 1.588125352, (-1)*0.39*x^2 + 1.459*x - y, 1.588125352 < x and x < 4, (x - 1.81)^2 + (y - 0.42)^2 + (-1)*0.94^2);

plots:-implicitplot(f, x = 0 .. 3, y = 0 .. 1.5, scaling = constrained,gridrefine=4,signchange=false);

In order to see that a sign change takes place you could try this:

plot3d([f,0],x = 0 .. 3, y = 0 .. 1.5);

Try this:

restart;
interface(typesetting=extended);
exp(-a); # No visible minus sign
interface(typesetting=standard);
exp(-a); #Fine also in Maple 2023.2

This produces the same as your code:
 

restart;
X := Matrix(
   [[1, -X3_3/2 - 1/2, 0, -X2_3], [-X3_3/2 - 1/2, -2*X3_4 - 1, X2_3, 0], [0, X2_3, X3_3, X3_4], [-X2_3, 0, X3_4, 1]]
            );
####
vars := [X3_4, X3_3, X2_3];
w := A^3 - A;

rts:=solve(w,A);
Pols := [(-A^2 + 1)/(3*A^2 - 1), (-A^2 - 1)/(3*A^2 - 1), A*(3*A^2 - 1)*1/(3*A^2 - 1)];
####
vals:=map2(eval,Pols,A=~{rts});
ecs:=[seq](vars=~(vals[i]),i=1..3);
Nelem := [seq(k, k = 1 .. numelems(vals))];
####
Xs:=[seq(eval(X, ecs[i]),i=Nelem)];
####
with(LinearAlgebra):

Eigs := Eigenvalues~(Xs);

BoolEigs:=map((y->is(Im(y) = 0 and 0 <= Re(evalf(y))))~ ,Eigs);

Maple can do it for you. It uses a finite difference method.
 

restart;
eq1 := (diff(f(x), x, x, x))*(a*beta*f(x)^2-1)+(diff(f(x), x))^2-2*a*beta*f(x)*(diff(f(x), x))*(diff(f(x), x, x))+(diff(f(x), x))*(M+k[1])-(diff(f(x), x, x))*f(x)-(alpha*theta(x)+delta*phi(x))/rho = 0;

eq2 := -(diff(theta(x), x, x))*K[SB]*(Df-(Rd+k[hnf]/k[bf])/Pr)+N[t]*K[SB]*(diff(theta(x), x))^2-N[b]*(diff(theta(x), x))*(diff(phi(x), x))-(diff(f(x), x))*(diff(theta(x), x))-lambda*theta(x)-mu*Ec*(M*(diff(f(x), x))^2+(diff(f(x), x, x))^2) = 0;

eq3 := diff(phi(x), x, x)+Le*Sr*(diff(theta(x), x, x))+Le*f(x)*(diff(phi(x), x)) = 0;

ics := f(0) = 0, (D(f))(0) = 0, theta(0) = 1, phi(0) = 1;

bcs := (D(f))(100) = 0, theta(100) = 0, phi(100) = 0;


Parameters1 := rho = 2063.905, k[hnf] = .29942, k[bf] = .2520, mu = .38694, a = .1, beta = 5, k[1] = 2.0, M = 10, alpha = 20, delta = 20, K[SB] = .5, Df = 3, Pr = 1.2, Rd = 5, N[t] = 1.2, N[b] = 1.0, lambda = 1.5, Ec = 5, Le = .1, Sr = .1;

 

sys:={eq1,eq2,eq3,ics,bcs};
res:=dsolve(eval(sys,{Parameters1}),numeric);
with(plots):
odeplot(res,[x,f(x)]);
scene:=[  [x,f(x)],[x,theta(x)], [x,phi(x)]  ];
display(Array([seq(odeplot(res,scene[i],thickness=2,color=blue),i=1..3)]));
scenediff:=evalindets(scene,function,fu->diff(fu,x));
display(Array([seq(odeplot(res,scenediff[i],thickness=2,color=red),i=1..3)]));

You may want a better look at the first part of the interval x = 0..100.
Just add the interval as done here for x = 0..20:
 

display(Array([seq(odeplot(res,scene[i],0..20,thickness=2,color=blue),i=1..3)]));
display(Array([seq(odeplot(res,scenediff[i],0..20,thickness=2,color=red),i=1..3)]));

restart;
Digits:=40;
plot( x*exp(-x),x=0..20,adaptive=true); 
?adaptive
# Quote: The adaptive = geometric option is used for non-parametric plots in the Cartesian coordinate system

The terrible result you get is caused by the default adaptive=geometric.

You could look at this:

restart;
showstat(`convert/rational`);
interface(verboseproc=2);
eval(`convert/rational`);