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

The following appears to work. Notice that abserr has been increased a lot.
 

A1 := dsolve(subs(para,rho=2, {bc, eq1, eq2}), numeric,method = bvp[midrich],
   initmesh=1024,maxmesh=12500, output=Array([seq( 0.01*i, i=0..100*N)]),abserr=1.5e-2):

You can use D or diff as you please in the ode (in fact you can also use the inert Diff or %diff:
General advice: Dont use N[0] when N is in use for something else. Here use N0 or N__0.
 

dsolve({D(N)(t)=r*N(t)-k*N(t)^2,N(0)=N0});
dsolve({diff(N(t),t)=r*N(t)-k*N(t)^2,N(0)=N0});
dsolve({Diff(N(t),t)=r*N(t)-k*N(t)^2,N(0)=N0});
dsolve({%diff(N(t),t)=r*N(t)-k*N(t)^2,N(0)=N0});

I have omitted the unknown functions name as the second variable. dsolve can handle it itself (generally).
Where you really need D is in giving initial conditions to higher order odes. See the help.

simplify has made use of the form simplify(expr, assume= ...)  before assuming was ever introduced.

They don't work quite the same:
 

restart;
f:=x->a*x;
simplify(sqrt(f(a)),assume=[a<0]);
simplify(sqrt(f(a))) assuming a<0;

The assume = [a<0] version works because simplify evaluates its arguments before it starts its general business as is the normal behavior in Maple.
The assuming version doesn't work here because it really works (as Carl points out) like this:
 

`assuming`([simplify(sqrt(f(a)))],[a<0]);

`assuming`  doesn't evaluate its first argument until the assumption a<0 is made on the visible a, i.e. the other a resulting from evaluating f(a) is not seen at that time, so internally simplify will then be looking at sqrt( a*a~), which it cannot do anything with since a and a~ are different at that time. At exit a~ is replaced by a, and you just have sqrt(a^2).

To use dsolve/numeric/bvp for your new problem you can use the symmetry.
 

restart;
t := -1:
S := 1:
R := 9:

EQ:={(diff(F(x), x $ 4)) + R*( diff(F(x),x)*diff(F(x),x$2)- F(x)*diff(F(x),x$3) ) + R*t*( F(x)*diff(F(x),x$5) - diff(F(x),x)*diff(F(x),x$4) ) - R*S^2*diff(F(x),x$2) =0};


IC:={D(F)(-1)=0, D(F)(1)=0,F(-1)=-1,F(1)=1,D(D(F))(1)=0};

(*Since it follows from EQ and your proposed boundary conditions that there is (could be) an odd solution F, i.e. one that satisfies
F(x) = -F(-x), you can restrict yourself to the interval -1..0. For an odd function F we have F(0) = 0, (D@@2)(F)(0)=0,(D@@4)(F)(0)=0, so we take: *)
 
bcs1:={D(F)(-1)=0,F(0)=0,F(-1)=-1,(D@@2)(F)(0)=0,(D@@4)(F)(0)=0};
sol1:=dsolve(EQ union bcs1,numeric,method=bvp[midrich]);
plots:-odeplot(sol1,[x,F(x)]);
A1:=Array([-1,-0.9,-0.8,-0.7,-0.6,-0.5,-0.4,-0.3,-0.2,-0.1,0]);
solA1:=dsolve(EQ union bcs1, numeric, method=bvp[midrich],output=A1);
M1:=solA1[2,1];
solA1[1,1]
M2:=Matrix(10,6,(i,j)->`if`(j=1,-M1[11-i,1],(-1)^(j+1)*M1[11-i,j]));
M:=< M1,M2>;
plot(M[..,1..2]);

You can also use odeplot, but then you need to make a structure like solA1:
 

sol:=Matrix(2,1);
sol[1,1]:=solA1[1,1];
sol[2,1]:=M;
plots:-odeplot(sol,[x,F(x)],-1..1);

You can do as follows. The names chosen by me can obviously be changed to your liking:
 

restart;
ode1:=diff(u(t), t) = a*u(t) - alpha__1 *u(t)*v(t)/(1 + k*u(t));
ode2:=diff(v(t), t) = -b*v(t) + alpha__1*u(t)*v(t)/(k*u(t) + 1) - alpha__2*v(t)*w(t);
ode3:=diff(w(t), t) = -c*(w(t) - wstar) + alpha__2*v(t)*w(t);
sys1:=PDEtools:-dchange({u(t)=U*p(tau),v(t)=V*q(tau),w(t)=W*r(tau),t=T*tau},{ode1,ode2,ode3},[p,q,r,tau]);
sys2:=solve(sys1,diff~({p,q,r}(tau),tau));
sys3:=expand(sys2);
# Now we pick T, U, and V.
# Looking at the first term on the right of the first ode in sys3 we see from the denominator that U must be chosen so k*U is dimensionless. 
# An obvious choice is U=1/k.
# After that T is chosen as 1/a. Then that first term is handled and is dimensionless.
# Moving to the next term in the same ode we choose W=a/alpha__2:
sys4:=eval(sys3,{U=1/k,T=1/a,V=a/alpha__1,W=a/alpha__2});
# Next step: Find the dimensionless groups and name them:
G:=[b/a,alpha__1/(a*k),c/a,alpha__2*c*wstar/a^2,alpha__2/alpha__1]=~[beta,delta,eta,zeta,a21];
SBS:=solve(G,{a,b,k,c,alpha__1});
## The dimensionless system with the dimensionless parameters beta,delta,eta,zeta,a21:
sys:=subs(SBS,sys4);

The system found is:
sys := {diff(p(tau), tau) = p(tau)^2/(1 + p(tau)) - p(tau)*q(tau)/(1 + p(tau)) + p(tau)/(1 + p(tau)), diff(q(tau), tau) = -q(tau)*r(tau)*p(tau)/(1 + p(tau)) - q(tau)*p(tau)*beta/(1 + p(tau)) - q(tau)*r(tau)/(1 + p(tau)) + q(tau)*p(tau)*delta/(1 + p(tau)) - q(tau)*beta/(1 + p(tau)), diff(r(tau), tau) = a21*q(tau)*r(tau) - eta*r(tau) + zeta}

Separate the text reion from the input region with F3.

Put the cursor in front of de:= ...  and press F3.
After that you could select the whole text line and convert to plain text.
In fact this is all you need to do.

Increase maxmesh. In this case maxmesh = 256 is sufficient.
 

sol:= dsolve(EQ union IC,numeric,maxmesh=256,output=Array([-1,-0.9,-0.8,-0.7,-0.6,-0.5,-0.4,-0.3,-0.2,-0.1,0,0.1,0.2,0.3,0.4,0.5,0.6,0.7,0.8,0.9,1]));

The sum to infinity can be made finite for plotting purposes:
 

restart;

pde := diff(T(x, y), x $ 2) + diff(T(x, y), y $ 2) = 0;
bc := T(0, y) = T1, T(a, y) = T2, T(x, 0) = T2, k*D[2](T)(x, b) = h*(-T(x, b) + T3);

sol1 := simplify((pdsolve([pde, bc], T(x, y)) assuming (0 < a, 0 < b)));

#and the constant value are:  a = 250, b = 4, k = 2.091, T1 = -5, T2 = 0, h = 100, T3 = 1000

#the plot range is 0<x<250, 0<y<4.
params:={a = 250, b = 4, k = 2.091, T1 = -5, T2 = 0, h = 100, T3 = 1000};
sol1n:=eval(sol1,params union {infinity=100}); # Choosing infinity =100 (!)

plot3d(rhs(sol1n),x=0..250,y=0..4);

### An animation in infinity (an attempt at justifying infinity = 100):
plots:-animate(plot3d,[eval(rhs(sol1),params),x=0..250,y=0..4],infinity=10..100, frames=9);

The animation in the parameter 'infinity':

Notice the difference between subs and eval:
 

subs(x=0,exp(x));  #exp(0)
%; # 1
eval(exp(x),x=0); # 1

That is by design.

Using assumptions helps:
 

restart;

expr:=a^2*(2*a^2*p^3-a^2*((p^2+1)^2*(a^2-1))^(1/2)-((p^2+1)^2*(a^2-1))^(1/2)*p^2+2*p*a^2-2*p^3+((p^2+1)^2*(a^2-1))^(1/2)-2*p)/((p^2+1)^2*(a^2-1))^(1/2)/(p^3-((p^2+1)^2*(a^2-1))^(1/2)+p)/(a^2-p^2-1);

expr2:=simplify(expr) assuming a>0,p::real;
res:=int(expr2,p);
diff(res,p)-expr2;
simplify(%) assuming a>0,p::real;
#####
idts:=indets(res,name);
select(type,idts,'indexed');
op~(0,%);
indets(%,`local`);
c:=op(%);
resC:=subs(c[5]=C,res);

You see that the result res contains an escaped local 'a' in an indexed form as a[5] besides your global a.

Had you used a different name e.g. b this escaped local would still have been 'a'..

Anyway, even this shouldn't happen.

Clearly a bug.

restart;
ode:=(y(x)-x*diff(y(x),x))/(y(x)^2+diff(y(x),x))=(y(x)-x*diff(y(x),x))/(1+x^2*diff(y(x),x));
sol:=dsolve(ode);
### Isolating y'(x):
odes:=solve(ode,{diff(y(x),x)});
dsolve~([odes]); # Fine