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

PS. I see that Carl Love beat me to it, but I shall leave my version. It doesn't hurt with several examples especially when you are new to Maple.

Could you not use events?
Here is a suggestion:
 

sysE := {
         diff(x(t), t) = v(t),
         diff(v(t), t) = t - (b(t)*f+(1-b(t))*(-f)),
         x(0) = 0,
         v(0) = 0,
         b(0) = 0
       };
solE := dsolve(sysE, numeric, parameters=[f],discrete_variables=[b(t)::boolean], events=[[v(t)=0,b(t)=1-b(t)]]);
solE(parameters=[1]):
solE(1);
odeplot(solE,[t,v(t)],0..100);

In the revised version of U it pays to reverse the order of integration:

restart;
#HAM := [1, 2, 3, 4, 5, 6]; 
U := q^6*sin(p); 
alpha := 1/3;
## Notice that the order of integration has been reversed:
for i from 1 to 6 do
   rho := i;
   f[i] := int(int(U/((x^rho-p^rho)^alpha*(y^rho-q^rho)^alpha), p = 0 .. x), q = 0 .. y) assuming x>0,y>0
end do;

Here as we saw we have D(theta)(0) = -1 (the second argument to Rproc):

This code:
 

plot3d(Rproc(thetaw,Dth0,initmesh=256,approxsoln=[f(eta)=0.5*tanh(eta),g(eta)=0.2*tanh(eta),theta(eta)=exp(-12*eta),R=1]),
     thetaw=1..2,Dth0=-3..-1,labels=[thetaw,Dth0,R]);

produces (after a while):


and this code produces an animation:
 

plots:-animate(plot,
     [Rproc(thetaw,Dth0,initmesh=256,approxsoln=[f(eta)=0.5*tanh(eta),g(eta)=0.2*tanh(eta),theta(eta)=exp(-12*eta),R=1]),
     thetaw=1..2,labels=[thetaw,R] ],
     Dth0=-3..-1,frames=10,trace=9);

Download MaplePrimes20-10-26_odebvp2.mw

Here is a start that may or may not with extra optional arguments to dsolve lead to something:
 

restart;
sys:=P*diff(f(eta), eta, eta, eta)/Q + lambda*theta(eta)*(S0 + S1*Qc*theta(eta)) + (n + 1)/2*(f(eta) + g(eta))*diff(f(eta), eta, eta) - n*diff(f(eta), eta)*(diff(f(eta), eta) + diff(g(eta), eta)) = 0, P*diff(g(eta), eta, eta, eta)/Q + (n + 1)/2*(f(eta) + g(eta))*diff(g(eta), eta, eta) - n*diff(g(eta), eta)*(diff(f(eta), eta) + diff(g(eta), eta)) = 0, (S2 - 2*R)*diff(theta(eta), eta, eta)/(S4*Pr) + (3*R)/(2*S4*Pr*(thetaw - 1))*diff((1 + (thetaw - 1)*theta(eta))^2, eta, eta) + (n + 1)/2*(f(eta) + g(eta))*diff(theta(eta), eta) = 0, f(0) = 0, g(0) = 0, theta(0) = 1, D(f)(0) = 1, D(g)(0) = c, D(f)(3) = 0, D(g)(3) = 0, theta(3) = 0;


sys[3]; # Term 2 has a removable singularity at thetaw=1
sys3:=map(normal,lhs(sys[3]))=0; # No singularity at thetaw=1
SYS:=sys[1..2],sys3,sys[4..-1];
P := 0.7952865489: Q := 1.239340734: S4 := 1.007143462:S0 := 0.8330195568:S1 := 0.8330195568: 
S2 := 1.087762421: Pr := 0.0157*2415/0.252: n:=3: lambda:=20: Qc:=2: c:=0.5: 
for s in SYS do s end do; # Viewing the system
Digits:=15:
#Boundary value problems of this kind are notoriously difficult, 
#so we begin with simply providing simple values for thetaw and R:
res:=dsolve(eval({SYS},{thetaw=1,R=2}),numeric,initmesh=256);
res(0);
dth0:=subs(res(0),diff(theta(eta),eta));
#The following procedure takes values of thetaw and D(theta)(0) as input and outputs the 
#corresponding value for R.
#Notice the extra boundary condition D(theta)(0)=Dth0:
Rproc:=proc(thetaw0,Dth0) local res;
   res:=dsolve(eval({SYS,D(theta)(0)=Dth0},thetaw=thetaw0),numeric,_rest);
   subs(res(0),R)
end proc;
#Here we should like to see the result 2 for R, but unfortunately we get an error instead. 
#Equally unfortunate is it that I'm not surprised:
Rproc(1,dth0,initmesh=1024);

Since you don't worry about the legitimacy of the simplification you can use the assumption  n::posint.
So you could do:
 

subsindets(expr,'exp(anything)^anything',f->simplify(f,exp)) assuming n::posint;

The following simple example is a reminder that the simplification is not always true. One should keep in mind that Maple uses the principal branch of the logarithm.
 

x:=2*Pi*I;
n:=-1+I;
exp(x)^n;        #  1
evalc(exp(x*n)); #  exp(-2*Pi)

dsolve doesn't use LinearAlgebra:-dAlambertianSolver, which is ised in DEtools:-particularsol.
That solver hangs. Why I don't know.
Try this:
 

restart;
ode:=diff(y(x),x$7)-2*diff(y(x),x$6)+9*diff(y(x),x$5)-16*diff(y(x),x$4)+24*diff(y(x),x$3)-32*diff(y(x),x$2)+16*diff(y(x),x)=exp(2*x)+x*sin(x)+x^2;
debug(LinearOperators:-dAlembertianSolver);
sol:=dsolve(ode); 
DEtools:-particularsol(ode,y(x));

Clearly a bug. But omitting the initial condition makes odetest work:
 

restart;
ode:=diff(y(x), x) - 2*y(x) = 2*sqrt(y(x));
ic:=y(0)=1;
maple_sol := 1 - 2*exp(x) + sqrt(y(x)) = 0;
my_sol:= y(x)^(1/2) = -1+2*exp(x);
odetest(maple_sol,ode) assuming x>0; # 0
odetest(my_sol,ode) assuming x>0 ;   # 0

This bug is very old. It is in Maple 8 too.
It appears that the bug is in `odetest/sysODE`, which is not used when ic is excluded.

Your problem just consists in finding an antiderivative P(x) of the right hand side satilfying the condition P(x) -> 0 for x -> -infinity:
 

restart;
sigma := x^2 + 1;
M:=1/10:
N:=1/2:
ode := diff(P(x), x) = -3*(M - 1)/(2*sigma^2) - 3*N/sigma^3;
sol:=int(rhs(ode),x=-infinity..x); # abuse of notation: x used in 2 ways
simplify(diff(sol,x)-rhs(ode)); # Check

Notice that I have used rationals instead of floating points.
The abuse of notation is allowed in Maple these days.

PS. Maybe you are not really interested in finding P(x), but are using this as training for odes having dependence on P(x) and not solvable by symbolic methods?