As far as producing [a,b,c,d,f''(0)] for various values of a,b,c,d you could do:

restart;
eq1:=(a,b)->{diff(f(eta),eta$3)-a*diff(f(eta),eta$1)^2+b*f(eta)*diff(f(eta),eta$2)=0};
bc:=(c,d)->{f(0)=0,D(f)(0)=1+c*(D@@2)(f)(0),D(f)(8)=d};
sol:=proc(a,b,c,d) [a,b,c,d,rhs(dsolve(eq1(a,b) union bc(c,d),numeric)(0)[-1])] end proc;
sol(1,2,1,2);

I ran your code with N=1000 and M=10 to save time. But my pdf file turned out the same as yours. 
Then I tried to Export/Encapsulated Postscript.
That looks fine and so does the pdf file I got from converting in ghostview/ghostscript.

I just tried changing the color to red before exporting to pdf. That worked!

Occasionally I have experienced the same in Maple 16 and 17. I use Standard Maple, worksheet mode, and Maple input (aka 1D input) almost exclusively.

In an existing worksheet after having executed the previous lines I get to a line (already written), the cursor is all the way to the left. Then it has happened (infrequently, yes) that Enter has no effect: Nothing happens. I always manage to get out of the problem somehow, using the mouse or whatever.
Since it happens infrequently and I have no idea of how to reproduce it I haven't mentioned it in this forum, nor have I filed an SCR (Software Change Request).
Surely it couldn't be a site license problem. My license is personal.

restart;
y1:=t->1/(4*cosh(t)^2);
I1:=int(y1(t)^2,t=-T/2..T/2);
convert(I1,exp);
res1:=combine(expand(%));
#Trying for a coarse approximation:
res:=simplify(eval(res1,{exp(-T/2)=0,exp(-3/2*T)=0}));
MultiSeries:-asympt(I1,T,10);
#Those agree.
#Now change variable:
eval(res1,T=2*ln(x));
asympt(%,x,10);
subs(x=exp(T/2),%);
combine(%);

You already got the answers to your problem, but I couldn't resist taking a look at your system.
I tried to change the parameter epsilon. The most convenient way to do that is by using the parameters option in dsolve/numeric:
restart;
sys:=epsilon*(diff(x(t), t)) = x(t)+y(t)-q*x(t)^2-x(t)*y(t),
                diff(y(t), t) = h*z(t)-y(t)-x(t)*y(t),
                p*(diff(z(t), t)) = x(t)-z(t);
#epsilon := 0.3e-1;#Commented out
h := .75;
p := 2;
q := 0.6e-2;
ics:=x(0) = 100, y(0) = 1, z(0) = 10;
sol := dsolve( {sys,ics}, type = numeric,maxfun=0,parameters=[epsilon]);
#Setting the parameter epsilon:
sol(parameters=[ 0.3e-1]); #The given epsilon
#An animation in phase space:
plots:-odeplot(sol,[x(t),y(t),z(t)],40..50,frames=30,numpoints=1000);
#Taking epsilon = 1
sol(parameters=[ 1]);
plots:-odeplot(sol,[x(t),y(t),z(t)],0..100,numpoints=1000);
#And try this animation in epsilon
p:=proc(eps) sol(parameters=[eps]); plots:-odeplot(sol,[x(t),y(t),z(t)],15..100,numpoints=1500) end proc;
plots:-animate(p,[epsilon],epsilon=.003..1,frames=100);

The solution is periodic, so you need only find the period.
Using Maple we can find that r(phi) satisfies
eq:=(6*r(phi)^5-2*r(phi)^3+r(phi)^2+(diff(r(phi), phi))^2)/r(phi)^4 = 13/4;
#Now let eq2 be
eq2:=subs(diff(r(phi),phi)=rp,r(phi)=r,eq);
#This is the equation of a closed curve, which is illustrated here:
plots:-implicitplot(eq2,r=0..1,rp=-3..3,gridrefine=3,crossingrefine=1);
#Using output=listprocedure
p:=dsolve({DE,ics},numeric,output=listprocedure);
R:=subs(p,r(phi));
PHI:=fsolve(R(phi)=.66666666,phi=5); #Using 2/3 exactly is problematic because of roundoff
#PHI is the period, as is illustrated here:
plots:-odeplot(p,[[phi,r(phi)],[phi,2/3]],0..PHI);
plots:-odeplot(p,[r(phi),diff(r(phi),phi)],0..PHI); #Phase plane plot

Using events agrees with the PHI we found:
p:=dsolve({DE,ics},numeric,output=listprocedure, events=[[[r(phi)-.6666666,diff(r(phi),phi)>0],halt]]);
p(5);

You have a system of pdes.
This should solve the system numerically:

restart;
PDE1:=diff(x(z,t),t)=(a/Pe)*diff(x(z,t),z$2)-a*diff(x(z,t),z)-(1-epsilon)/epsilon*3*Bi*(x(z,t)-1)/(1-Bi*(1-1/ksi(z,t)));
a:=2:Pe:=3:Bi:=5:epsilon:=0.85: ##Assignment
PDE2 := diff(ksi(z,t),t) = (b*Bi*x(z,t)-1)/ksi(z,t)^2/(1-Bi*(1-1/ksi(z,t)));
IC1:=x(z,0)=0; ## c = x?
IC2:=ksi(z,0)=1;
bc2:=D[1](x)(1,t)=0;
bc1:=x(0,t)-1/Pe*D[1](x)(0,t)=0;
b:=1:  ###Not given
sol:=pdsolve({PDE1,PDE2},{IC1,IC2,bc1,bc2},numeric);
sol:-plot([x(z,t),[ksi(z,t),color=blue]],t=.3);
sol:-animate([x(z,t),[ksi(z,t),color=blue]],t=.3);

Here are two ways.

1. Make separate plots. Then display them together.
p1:=plot(r^2,r=0.2..0.5):
p2:=plot(sin(r),r=0.5..2):
plots:-display(p1,p2);

2. You could use piecewise since the intervals don't overlap.
plot(piecewise(r<0.5,r^2,sin(r)),r=0.2..2,discont=true);

You could split as you do by hand:

dsolve({eq3=0,phi(0)=0,phi(h)=1},phi(eta));
res1:=eval(%,{bcs});
eq:=simplify(eval(eq2=0,res1));
res_T:=dsolve({eq,T(0)=0,T(h)=1});
res_phi:=eval(res1,res_T);
odetest([res_T,res_phi],[sys_ode]);

The result you show contains an unevaluated integral. In order for dsolve to do its job it had to do some integration. _U1 is just an integration variable needed in the integral. The underscore indicates that dsolve had to come up with the name itself. It is a global variable, so if you like you could replace it with some other name (say xi) by doing subs(_U1=xi, o[cx](t).
Anyway, Maple didn't come up with a finished result.
To get any further we need to have the Maple worksheet.

With the correct syntax as Markiyan is pointing out there is no problem for dsolve.
However, the output is extremely long.
method=laplace gives a much shorter output:

restart;
sys_ode:=diff(A(t), t) = -k1*A(t)+k2*B(t),
               diff(B(t), t) = k1*A(t)+(-k2-k3)*B(t)+k4*C(t),
               diff(C(t), t) = k3*B(t)-k4*C(t);
ics := A(0) = S, B(0) = 0, C(0) = 0;
res:=dsolve({ics, sys_ode},[A(t),B(t),C(t)],method=laplace);
#We can make it look shorter if we use
indets(res,radical);
op~(1,%);
R:=op(%);
subs(R=rr,res);

For general rate constants k1, ... ,k4 the result is messy, however.

Since you are dealing with foxes and rabbits I suspect that you have a sign error in your equation for dr/dt.
The rabbits get eaten by the foxes (at least killed) so since alpha > 0 you must have
diff(r, t) = 2*r - alpha*r*f;

Assuming this then you have a classical Volterra-Lotka system, where all solutions are periodic.
Now try:
restart;
DE := diff(r(t), t) = 2*r(t)-alpha*r(t)*f(t), diff(f(t), t) = -f(t)+alpha*r(t)*f(t);
params := alpha = .3;
initv := r(0) = 101, f(0) = 2;
EQ := [op(subs(params, [DE])), initv];
 
EQ1 := dsolve(EQ, numeric);

plots:-odeplot(EQ1,[r(t),f(t)], t = 0 .. 100, numpoints=5000,axes = frame, color = green, title = "Rabbit and Fox") ;
#Notice the spiraling behavior: That should not be there! All solutions are periodic.
The problem is numerical. Small errors unfortunately build up. If you trust that the solution is periodic then just don't use that large a time interval, 0..20 should be sufficient.

It requires a little work to see that solutions are periodic, but here is an indication and illustration.

Eliminating t we get this ode for f(r):
eq:=diff(f(r),r)=(-f(r)+alpha*r*f(r))/(2*r-alpha*r*f(r));
eq1:=dsolve(eq,implicit);
eval(eq1,{r=101,f(r)=2,params});
solve(%,{_C1});
eq2:=eval(eq1,% union {params,f(r)=f});
sol:=solve(eq2,f);
#Actually the first of the two solutions found is wrong, but the second is correct:
sol1:=sol[2];
#To see that first is wrong do
eval(eq2,f=sol[1]);
simplify(%);
eval(%,r=10.); # This ought to have been 0=0 at least approximately.
#The correct other soltion uses another branch of LambertW:
sol2:=eval(sol[2],LambertW=curry(LambertW,-1));
plot([sol1,sol2],r=0..105);

If I understand you correctly, then x11 lists the x1-values corresponding to the times listed in 'tim'. And similarly with y11 and z11.
Using Optimization:-Minimize instead of NonlinearFit:

sol:=dsolve({a1,b1,c1,ICS}, numeric, method=rkf45, parameters=[k1,k2,k3,k4,k5,k6,k7,k8,k9],output=listprocedure);
X,Y,Z:=op(subs(sol,[x1(t),y1(t),z1(t)]));
tim := [seq(n, n=1..27)];
N:=nops(tim):
ans:=proc(k1,k2,k3,k4,k5,k6,k7,k8,k9) sol(parameters=[k1,k2,k3,k4,k5,k6,k7,k8,k9]);
 add((X(tim[i])-x11[i])^2,i=1..N)+add((Y(tim[i])-y11[i])^2,i=1..N)+add((Z(tim[i])-z11[i])^2,i=1..N)
 end proc;
ans(.001,.002,.003,.001,.002,.003,.001,.002,.003);
Optimization:-Minimize(ans,initialpoint=[.001,.002,.003,.001,.002,.003,.001,.002,.003]);
#The error reported (i.e. the sum of squares of residuals) is with the given initialpoint 6.5*10^12. That is huge.
I doubt you can do substantially better with the data you got. Do you think that the linear model is any good?

1. You forgot to assign to a.
2. Use unapply in defining psi_5 .. psi_8:
psi_5:=unapply(psi_1(x,y)*chi(x,y),x,y);
psi_6:=unapply(psi_2(x,y)*chi(x,y),x,y);
psi_7:=unapply(psi_3(x,y)*chi(x,y),x,y);
psi_8:=unapply(psi_4(x,y)*chi(x,y),x,y);

3. Wrap the definitions of the relevant matrices in evalf so that the integrals get computed.
To see if E is a numerical matrix you can do:
interface(rtablesize=infinity);
E;

4. The computation should be fast now, but notice that the datatype of E is anything:
op(E);
To change that you could do E:=Matrix(E,datatype=float);

The results of doing what I described in my last comment above made me think that an approximate solution for the system in the scaled variables y and psi would have
psi[1](t) = N-t
psi[3](t) = 0
y[1](t) = 1
y[2](t) = 1

with psi[2] and y[3] to be determined.
Inserting the 4 "known" functions into the system one finds very simple differential equations for psi[2] and y[3].
Using the conditions psi[2](N) = 0 and y[3](0)=1 dsolve (without 'numeric') finds easily that
psi[3](t) = 0
and
y[3](t) = 2.1891-1.1891*exp(-.67*t).

OK, so how good is this approximate solution?
As far as satisfying the system of odes measured by the difference between rhs and lhs of the equations the answer is that the largest difference when N=100 is 0.0018.

MaplePrimes13-12-07b.mw