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

Two exact solutions are easily found:

restart;
###
A1 := 8*Pi^3*R^2*n(x)^4*m+(2*Pi*sin((1/2)*x)*m*omega0*p+Pi*sin((1/2)*x)*m*omega0+3*Pi^2*(diff(n(x), x, x)))*n(x)^3+(-2*sin((1/2)*x)^2*m^2*omega0^2*p^2+2*cos((1/2)*x)^2*m^2*omega0^2*p^2-2*sin((1/2)*x)^2*m^2*omega0^2*p+2*cos((1/2)*x)^2*m^2*omega0^2*p)*n(x)^2+(-4*(diff(n(x), x$2))*sin((1/2)*x)^2*m^2*omega0^2*p^2-8*sin((1/2)*x)*(diff(n(x), x))*cos((1/2)*x)*m^2*omega0^2*p^2-4*(diff(n(x), x$2))*sin((1/2)*x)^2*m^2*omega0^2*p-8*sin((1/2)*x)*(diff(n(x), x))*cos((1/2)*x)*m^2*omega0^2*p)*n(x)+8*sin((1/2)*x)^2*(diff(n(x), x))^2*m^2*omega0^2*p^2+8*sin((1/2)*x)^2*(diff(n(x), x))^2*m^2*omega0^2*p;

###
R := 1; m := 1; p := 10; omega0 := 1000;
###
bc:=n(0) = n(4*Pi), D(n)(0) = D(n)(4*Pi);
###
### Ansatz: n(x) = A*sin(x/2)
eval(A1,n(x)=A*sin(x/2));
factor(%);
Asol:=solve(%/A^2,A);
### Thus we have the two exact solutions:
sol1,sol2:=Asol[1]*sin(x/2),Asol[2]*sin(x/2);
odetest(n(x)=sol1,A1); # 0
odetest(n(x)=sol2,A1); # 0
### bc is clearly satisfied for both solutions. 
### Plot:
plot([sol1,sol2],x=0..4*Pi);

You can rather easily get a larger image by using the following code, where the appropriate size depends on your monitor's screensize:

DGR:=DrawGraph(GH, stylesheet = "legacy"):
plots:-display(DGR,scaling=unconstrained,size=[1800,default]);

You can also replace 'default' by an appropriate number (e.g. 800).

You should certainly also look into the many options for DrawGraph.

The following copy has luckily been shrunk by MaplePrimes:

You can do like this:

restart;
Digits:=15:
sigma := 1 + x^2/2;
k := 0.1;
Q := 0.5516;
### Use unapply:
lambda1 := unapply(-3*(7*k^3*sigma^3 + 32*Q*k^2*sigma^2 - 11*k^2*sigma^3 + 54*Q^2*k*sigma - 44*Q*k*sigma^2 + 11*k*sigma^3 + 36*Q^3 - 54*Q^2*sigma + 32*Q*sigma^2 - 7*sigma^3)/(20*sigma^4),x);
###
N:=6/10/10^(-6);
M := Matrix(N + 1, 2, datatype = float[8]);
###
for i from 1 to N+1 do x:=0.6*(i-1)/N; M[i]:= <lambda1(x)| x > end do:
M;
M[-1]; # the last row

But this is much faster:
 

restart;
Digits:=15:
sigma := 1 + x^2/2;
k := 0.1;
Q := 0.5516;
lambda1 := unapply(-3*(7*k^3*sigma^3 + 32*Q*k^2*sigma^2 - 11*k^2*sigma^3 + 54*Q^2*k*sigma - 44*Q*k*sigma^2 + 11*k*sigma^3 + 36*Q^3 - 54*Q^2*sigma + 32*Q*sigma^2 - 7*sigma^3)/(20*sigma^4),x);
N:=6/10/10^(-6);
V:=Vector(N+1,i-> (i-1)*0.6/N,datatype=float[8]);
W:=evalhf(map(lambda1,V));
M:=<W|V>;

The following is simpler and works in Maple 12 too:
 

p:=piecewise(`ϑ`(tau)>=ϑl,1,0);
eq := diff(`ϑ`(tau), tau, tau)+6.666666666*sin(`ϑ`(tau))+66.66666666*cos(`ϑ`(tau))^2*p*(`ϑ`(tau)-.7227342478)+66.66666666*sin(`ϑ`(tau))^2*p*(`ϑ`(tau)-.7227342478);
`ϑl` := .7227342478;
ic := `ϑ`(0) = 0, (D(`ϑ`))(0) = 8;
ld := dsolve([eq, ic], numeric, range = 0 .. 5);
odeplot(ld,[tau,p],refine=1,thickness=2, color=blue);