In a vain attempt to speed up things I tried turning the bvp into an initial value problem with parameters. That would involve fsolve so the computation overall took about twice as much time.

However, in the proces I discovered that the boundary value problem can have at least two solutions for given values of a,b,c,d.This is certainly the case for a=b=c=d=0.1 as the following shows. The two solutions have different values for f''(0).

restart;
eq1:= diff(f(eta),eta$3) - a*diff(f(eta),eta$1)^2 + b*f(eta)*diff(f(eta),eta$2) = 0:
bc:= f(0)=0, D(f)(0) = 1+c*(D@@2)(f)(0), D(f)(8)=d:
sys:= unapply({eq1,bc}, a, b, c, d):
Sol:= (a,b,c,d)-> dsolve(sys(a,b,c,d), numeric):
#Example a = b = c = d = 0.1.
plots:-odeplot(Sol(.1,.1,.1,.1),[eta,f(eta)],0..8,caption="Solution 1"); p1a:=%:
#############
ic:= f(0)=0, D(f)(0) = 1+c*f2, (D@@2)(f)(0)=f2:
solpar:=dsolve({eq1,ic}, numeric,parameters=[a,b,c,f2],output=listprocedure):
F1:=subs(solpar,diff(f(eta),eta));
#We make a procedure returning f'(8)-d:
p:=proc(a,b,c,d,f2)
   if not type([_passed],list(numeric)) then return 'procname(_passed)' end if;
   solpar(parameters=[a,b,c,f2]);
   F1(8)-d
end proc;
plot(p(.1,.1,.1,.1,f2),f2=-1..1); #Two zeros!!
r1:=fsolve(p(.1,.1,.1,.1,f2)=0,f2=-.3);
r2:=fsolve(p(.1,.1,.1,.1,f2)=0,f2=-.8);
solpar(parameters=[.1,.1,.1,r1]):
plots:-odeplot(solpar,[eta,f(eta)],0..8,color=red,caption="Solution 1 again"); p1b:=%:
solpar(parameters=[.1,.1,.1,r2]):
plots:-odeplot(solpar,[eta,f(eta)],0..8,color=blue,caption="Solution 2"); p2b:=%:
#The second solution can also be produced by solving the bvp problem directly when an approximate solution is known, and of course we have one.

Sol2:= dsolve(sys(.1,.1,.1,.1), numeric,approxsoln=solpar);

plots:-odeplot(Sol2,[eta,f(eta)],0..8,caption="Solution 2 again"); p2a:=%:
plots:-display(p1a,p2a);
plots:-display(p1b,p2b);
plots:-display(p1a,p2a,p1b,p2b);

I don't quite understand what you mean by plotting epsilon against time.
But if the intention is to plot the dependence on epsilon of the solution y then you could do something like this:

restart;
dsys1 :=diff(y(t),t)=y(t)*((1-y(t)/3-epsilon)-0.8*y(t)/(y(t)^2+0.5^2));
ini1:=y(0)=0.5:
dsol1 :=dsolve({dsys1,ini1},numeric, output=listprocedure,range=0..1,parameters=[epsilon]);
dsolu:=subs(dsol1,y(t));
#Setting the parameter epsilon for a test plot:
dsol1(parameters=[.02]);
plots:-odeplot(dsolu,[t,y(t)],0..1,caption=typeset(epsilon = 0.02));
#Now set up a procedure which returns unevaluated for non-numeric input, but otherwise sets the parameter epsilon to the input value (eps) and simply returns the procedure dsolu.
p:=proc(eps)
  if not type(eps,numeric) then return 'procname(_passed)' end if;
  dsol1(parameters=[eps]);
  dsolu
end proc;
#An animation:
plots:-animate(plot,[p(epsilon)(t),t=0..1],epsilon=0..1);
#A 3d plot:
plot3d(p(epsilon)(t),t=0..1,epsilon=0..1);

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?