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?

I don't think that you can prove anything in Maple.
But try this:
 

pdsolve(diff(g(y,z),y)=0);

@acer You can simply replace all occrurrences of subs in `plots/animate` by eval with this simple code:

`plots/animate`:=subs(subs=((x,y)->eval(y,x)),eval(`plots/animate`)):

This is quite harmless and works. I have been doing it for years. The problem is present in versions up to and including Maple 14.
It is unnecessary in later version: they don't use subs, so even here it is harmless.
After that redefinition the following code works quite nicely:
 

animate (pointplot3d, [ [b(t),c(t),d(t)], symbol=box, color=blue],t=15..100,
background = odeplot(sol, [x(t),y(t),z(t)],t=15..100,numpoints=7000),frames=500,axes=boxed);

I just tried in Maple 12.

Here is a straightforward way:
 

restart;
p:=<1,2,3,4,5,6>;
vs:=<.1,.6,NA,-1.4,NA,6.7>;
p1:=Vector():
vs1:=Vector():
j:=1:
for i from 1 to numelems(p) do
  if vs[i]<>NA then p1(j):=p[i]; vs1(j):=vs[i]; j:=j+1 end if
end do;  
vs1;
p1;

There is indeed a singularity of your ode. No way of avoiding that unless your numbers are changed:

Download MaplePrimes20-08-26_ode_ic_sing.mw