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:

It seems that algsubs won't touch attempts to replace functions like int, Int, sin, or f in f(x):
 

restart;
expr:=int(f(x),x);
algsubs(f(x)=g(x),expr); # Works
algsubs(f=g,expr); # Doesn't do anything
algsubs(sin=cos,sin(x)); # Doesn't do anything
algsubs(x=y,expr); # Works

That problem with y(0) = 0 has infinitely many solutions (a well known fact):
 

restart;
ode := diff(y(x),x) = 2*sqrt(y(x));
p:=seq(piecewise(x<n,0,(x-n)^2),n=0..4);
plot([p],x=-1..5);
seq(odetest(y(x)=q,ode),q=p);

With 17 equations we may assume that your solution is numerical.

I give here a very simple example: Just one ode and it can be solved symbolically as well as numerically:
 

restart;
ode1:=diff(x(t),t)=-x(t)+sin(t);
ic1:=x(0)=1;
res1:=dsolve({ode1,ic1},numeric); # The numerical solution for x
plots:-odeplot(res1,[t,x(t)],0..20); 
### Now introduce the integral of x which is called y: 
INT:=y(t)=int(x(s),s=0..t);
### Differentiate that to get an extra ode:
ode2:=diff(INT,t);
ic2:=eval(INT,t=0); # The initial value for y
res2:=dsolve({ode1,ode2,ic1,ic2},numeric); # Numerical solutions for x and y
plots:-odeplot(res2,[[t,x(t)],[t,y(t)]],0..20); 
### Now for comparison we solve exactly (symbolically):
sol1:=dsolve({ode1,ic1});
map(int,sol1,t=0..T);
sol2:=y(t)=subs(T=t,rhs(%));
plot([rhs(sol1),rhs(sol2)],t=0..20); 
### The difference between the numerical and the exact y(t):
plots:-odeplot(res2,[t,y(t)-rhs(sol2)],0..20);