@AHSAN In your worksheet solve gives its answer as a RootOf expression.
That expression is not at all useless.
 

restart;
pg1 := -(888888889*sigma^6*(-(2*k*sigma + 3*Q - sigma)/sigma^2)^(5/2) - 888888889*sigma^6*((k*sigma + 3*Q - 2*sigma)/sigma^2)^(5/2) + 3333333333*(-(2*k*sigma + 3*Q - sigma)/sigma^2)^(3/2)*k*sigma^5 + 3333333333*((k*sigma + 3*Q - 2*sigma)/sigma^2)^(3/2)*k*sigma^5 + 6666666666*(-(2*k*sigma + 3*Q - sigma)/sigma^2)^(3/2)*Q*sigma^4 - 3333333333*(-(2*k*sigma + 3*Q - sigma)/sigma^2)^(3/2)*sigma^5 + 6666666666*((k*sigma + 3*Q - 2*sigma)/sigma^2)^(3/2)*Q*sigma^4 - 3333333333*((k*sigma + 3*Q - 2*sigma)/sigma^2)^(3/2)*sigma^5 + 15000000000*k^2*lambda*sigma^2 + 60000000000*Q*k*lambda*sigma - 30000000000*k*lambda*sigma^2 + 60000000000*Q^2*lambda - 60000000000*Q*lambda*sigma + 15000000000*lambda*sigma^2)/(5000000000*(k*sigma + 2*Q - sigma)^2*sigma^3);
S := solve(pg1 = 0, sigma);
Sa:=allvalues(S):
numelems({Sa});
evalf(eval([Sa],{k = 2, Q = 3, lambda = 4}));
select(type,%%,realcons);

While evalc assumes that any variables present are real, convert/trig doesn't.
In fact the equation exp(I*z) = cos(z) + I*sin(z) holds for all z in the complex plane.

evalc(exp(I*x));
convert(exp(I*x),trig);

Here you will see the solution method:
 

restart;

ode:=diff(y(x),x)*tan(diff(y(x),x))+ln(cos(diff(y(x),x))) = y(x);
infolevel[dsolve]:=5:
dsolve(ode);
DEtools:-odeadvisor(ode); # Wrong

The printout is:

Methods for first order ODEs:
-> Solving 1st order ODE of high degree, 1st attempt
trying 1st order WeierstrassP solution for high degree ODE
trying 1st order WeierstrassPPrime solution for high degree ODE
trying 1st order JacobiSN solution for high degree ODE
trying 1st order ODE linearizable_by_differentiation
trying differential order: 1; missing variables
<- differential order: 1; missing "x" successful

If you try

indets(sol,name);

you get {O,O, a, b, x, y}.
If you continue with
 

indets(%,`local`);

you get  {O,O}

so the two O's are different locals.
Three O's appear in the expression, but two of those are equal.
If you try this:

oh:=op(indets(denom(op(1,rhs(sol))),`local`));
subs(oh=XX,sol);

you will see which two ones are equal.

The complex values must come from eq3, where you have (1 + delta*theta(eta))^n and where later n = 0.2.
Just try

(-2)^0.2; # answer: 0.9293164906 + 0.6751879524*I (the principal 5th root.).
You could then think of using surd (see the help). That turns out not to work in dsolve/numeric/bvp because the round disappears and you get surd( ... , 5. ), i.e a float instead of an integer as the second argument.
You can use eq3_new instead of eq3:

eq3_new:=subs((1 + delta*theta(eta))^n=abs(1 + delta*theta(eta))^n*signum(1 + delta*theta(eta)),eq3);

The problem is illustrated here:

surd(-2, round(1/0.2)); #Here it works because the integer 5 remains an integer
evalf(%);
### The alternative
abs(-2)^0.2*signum(-2);

You run into singularity problems though, so there are obviously other problems with your system.
With these extra arguments to dsolve
method=bvp[midrich],initmesh=1024
you get the error:

Error, (in dsolve/numeric/bvp) matrix is singular

The exp term in eq3 gives problems here:
 

exp(-E/(1 + delta*theta(eta)));
####
plot(exp(-10/(1 + 10*theta)),theta=-0.5..0.5,thickness=3);

This is really just an observation, but it appears that the mere presence of sqrt(-x) or sqrt(-y) makes the difference:

restart;
A:=(x+y-2*sqrt(x*y))+sqrt(-x);
simplify(A) assuming negative;
restart;
A:=(x+y-2*sqrt(x*y))+sqrt(-x);
simplify(A) assuming positive;

Thus this trick:

restart;
# Example 2

B:=x+y-2*sqrt(x*y);

simplify([B]+sqrt(-x)) assuming negative;

I notice that solving for diff(y(x),x) and applying dsolve works (fast):
 

ode1,ode2:=solve(x^3*diff(y(x),x)^2+x*diff(y(x),x)-y(x) = 0,{diff(y(x),x)});
dsolve(ode1);
dsolve(ode2);

restart;
ode1:=diff(u(x), x, x) = u(x)*sin(x);
ode2:=diff(ode1,x);
bcs:=u(0)=2,u(1)=1,(D@@2)(u)(1)=sin(1);
Digits:=10: # The default
res:=dsolve({ode2,bcs},numeric,abserr=1e-13);
plots:-odeplot(res,[x,diff(u(x),x,x)-u(x)*sin(x)],0..1,caption=typeset(diff(u(x),x,x)-u(x)*sin(x)));

A comment or rather a rhetorical question:

If you were to design an algorithm for dsolve/numeric (bvp or ivp) that returned the highest derivative too, what would you make it return other than the right hand side?
#####
Looking at what `dsolve/numeric` does to an ivp like this:
 

restart;
ode1:=diff(u(x), x, x) = u(x)*sin(x);
ode2:=v(x)=diff(u(x),x,x);
ics:=u(0)=2,D(u)(0)=0;
#debug(`dsolve/numeric`);
res:=dsolve({ode1,ode2,ics},numeric,abserr=1e-13,relerr=1e-10);
plots:-odeplot(res,[x,v(x)-u(x)*sin(x)],0..1,caption=typeset(v(x)-u(x)*sin(x)));

If you remove the # then the debugging output reveals what dsolve/numeric does here:
 

INFO["proc"] := proc(N, X, Y, YP) local Z1; Z1 := Y[1]*sin(X); Y[3] := Z1; if N < 1 then return 0; end if; YP[2] := Z1; YP[1] := Y[2]; 0; end proc;

As you see the right hand side is kept in the local Z1 and Y[3] (i.e. v) gets that value as does YP[2] i.e. u''.
This simply means that the right and side is returned for v.

The reason sol is not changed by writing it in another order is that once in memory the order is taken from there.
Maple doesn't keep several versions of expressions that it recognizes as the very same. I think simply for the sake of efficiency.
Just try this:

restart;
b+a;
a+b;
b*a;
a*b;
restart;
a+b;
b+a;
a*b;
b*a;

The following has nothing to do with a proof, so if that is what you really want then it is worthless.
It does show, however, how in Maple you can come from sin(2*x) to 2*sin(x)*cos(x) and vice versa:

expand(sin(2*x));
combine(%);
                           

Assuming that we accept as known that exp(z+w) = exp(z)*exp(w) for all complex z and w, then you can begin to call the following a proof:
 

E2:=exp(2*I*x);
evalc(E2);
Im(%) assuming x::real;
expand(E2);
evalc(%);
Im(%) assuming x::real;

I just ran your worksheet beginning of course with restart.
I used Maple 2020. No problem!
In Maple 18 I had your problem. Although the local declaration of gamma in Maple 18 is valid, dsolve,numeric,parameters doesn't like it being local.
Solution: use another name than gamma.
pi doesn't have to be declared local, Pi would.

The problem exists also in Maple 2015 and 2016, but not in Maple 2017 and later.
###########################
Try replacing gamma with gama (or some other name).
Omit declaring gama local. (There is no reason to do that anyway).
Your worksheet runs.
Now try the following:
 

ans := dsolve( { ODE1, ODE2, ODE3, ODE4, ODE5, ODE6, ODE7,
                   B(0) = B0, C(0) = C0, E(0) = E0, G(0) = G0, H(0) = H0, J(0) = J0, K(0) = K0
               
                 },
                 parameters = params,
                 numeric
               );
ans(parameters);
type(gama,`local`); # false
type(B0,`local`); #false
lhs~(ans(parameters));
type~(%,`local`); # all true

This is the same in Maple 2020. I suppose the reason simply is that if you after having defined ans actually assign to any of the parameters then ans is not ruined.

You can do something like this:
 

EVTS:=[[t = 2000, [b(t) = 0, Q[1](t) = 0]], [t = 2500, b(t) = 1]];
initialConditions := Q[1](0) = 0., H[1](0) = 1.5, H[2](0) = 1.5, b(0) = 1;
momentumBalance[1] := diff(Q[1](t), t) = (pipefactor*(H[1](t) - H[2](t)) + pumphead(t) - frictionfactor*abs(Q[1](t))*Q[1](t))*b(t);

####
res := dsolve({initialConditions, height[1], height[2], momentumBalance[1]}, numeric, output = listprocedure, discrete_variables = [b(t)::boolean], events = EVTS);
####
plots:-odeplot(res,[t,Q[1](t)],0..nts);
plots:-odeplot(res,[[t,H[1](t)],[t,H[2](t)]],0..nts);

The plots:

A:=int(diff(u(x,t),x),x=-infinity..infinity,continuous);

Answer:  A := -Limit(u(_X, t), _X = -infinity) + Limit(u(_X, t), _X = infinity)

Thus by your assumption A = 0.

restart;
f:= 2/(3-x);
s:=convert(f,FormalPowerSeries);
op(s);
simplify(op(1,s)); 
eval(%,x=1);

You could use odeplot from the plots package.

Here is an example which I have put together hastily (I admit):
 

restart;
sys := {diff(x(t), t) = x(t) - y(t), diff(y(t), t) = 2*x(t) - y(t) - x(t)^2};
DEtools:-DEplot(sys,[x(t),y(t)],t=-2..5,x=-1..3,y=-2..3,[[x(0)=2,y(0)=3],[x(0)=1/2,y(0)=1/2]],linecolor=blue);
ics:={x(0)=2,y(0)=3};
res:=dsolve(sys union ics,numeric);
plots:-odeplot(res,[t,sqrt(x(t)^2+y(t)^2)],0..2);
res(initial);
res(initial=[0,1/2,1]);
plots:-odeplot(res,[t,sqrt(x(t)^2+y(t)^2)],0..8);
?dsolve,numeric,interactive
PL:=proc(IC::listlist,{scene::list:=[x(t),y(t)],range::range:=0..2}) local ic,p;
   for ic in IC do
      res('initial'=ic);
      p[ic]:=plots:-odeplot(res,scene,range,_rest)
   end do;
   plots:-display(seq(p[ic],ic in IC))
end proc;
###Examples:
PL([seq([0,i,1/2],i=-1..1,0.1)]);
PL([seq([0,i,1/2],i=-1..1,0.1)],range=-2.3..2);
PL([seq([0,i,1/2],i=-1..1,0.1)],range=-3..4,view=[-1..3,-2..3]);
PL([seq([0,i,1/2],i=-1..1,0.1)],range=-3..4,scene=[t,sqrt(x(t)^2+y(t)^2)],view=[-3..4,0..10]);

The last plot: