@Adri van der Meer You may want to accompany that by

M:=simplify(M) assuming real;

It would help a lot in getting people interested in your problem if they didn't have to enter the numbers themselves. You only provided an image.

@J4James If 8 is really a replacement for infinity, then one may speculate that the non-monotone solution of the ivp converges toward the monote solution on any bounded interval eta=0..E as "8" -> infinity.
If that is the case then the higher of the two f''(0) values would be the one to choose.

@Carl Love If so yes there are 2 solutions.
Try with 8 replaced by 12:
plot(p(.1,.1,.1,.1,f2),f2=-1..1,-.1..0.1); #Two zeros
r1:=fsolve(p(.1,.1,.1,.1,f2)=0,f2=-.3);
r2:=fsolve(p(.1,.1,.1,.1,f2)=0,f2=-.4);

@Carl Love Initially I just wanted to make use of the parameters option to see if that would save time. For that I had to turn the bvp into an ivp. Thus I needed f(0), f'(0), f''(0) given either concretely or as parameters.
The condition at eta=8 then had to be imposed afterwards.

@adel-00 I don't think we quite understand each other. Plotting beta versus t? Which t? I got the impression that that t was to be determined by the requirement (in this new case) N'(t) =0. 
What I was saying was that if for some t = t0 we had N'(t0) =0, then the right hand side of the ode (dsys1) would be zero, so either N(t0)=0 or N(t0) = a solution of a cubic (as in the previous problem).
So say we have N(t0) = N1 where N(t) = N1 makes the right hand side of the ode 0.
Then the ivp consisting of dsys1 and N(0)=N1 would have the constant solution N(t)=N1.
However, it also has the solution with initial condition N(0)=0.61299467707789. By uniqueness the solutions are the same, so we have that N1 =0.61299467707789. This determines beta. In this case we find:
subs(N(t)=rhs(ini1),rhs(dsys1));
solve(%=0,beta);
#returns beta roughly equal to zero.
My point, though, is that the problem as I understand it makes no sense.

@adel-00 You have a first order autonomous differential equation satisfying the conditions for uniqueness of the initial value problem. This implies that if y'(t) ever becomes 0 then it will have been zero at all times, i.e. that y is in fact constant.
For y' to be zero the right hand side of dsys1 has to be zero. The solutions to rhs = 0 are
y = 0 and the roots of the third order polynomial
pol:=-(1/3)*y^3+(1-epsilon)*y^2-(53/60)*y+(1-epsilon)/4
For some values of epsilon this polynomial has 3 real roots for other values only one.
You have y(0)=1/2 so if you want your solution to be constant that constant is 1/2.
Inserting y=1/2 into pol=0, you find that epsilon = 1/30.

@adel-00 This is an animation with epsilon as animation parameter. Each frame is a plot of y versus t.

This happens rather rarely and, as Markiyan points out, chances of getting an answer are greatly diminished. So I don't see any point in disallowing other languages.
However, I must admit to some irritation when somebody posts in a language other than English, in particular when it is in a language that I don't understand.

It is somewhat strange, yes, and maybe impolite to post a question in a language other than English on MaplePrimes. If the reader doesn't understand the language he is left out and cannot answer.
It doesn't happen very often though and, as Markiyan says, such a question is not likely to get many answers, if any at all. So I don't see the need for disallowing other languages.

@ozlem I am rather convinced that the system has a singularity at t = 0.31292.

Try with timestep=1e-5 this:
p:=sol:-value([x(z,t),ksi(z,t)]);
res:=p(1,.31292);
eval(rhs(PDE2),res); #returns 277.5303397
eval(rhs(PDE1),res);  # returns 241.5830864

@ozlem I cannot see the file you tried to udload, but the answer is (for t=0..0.3):

sol:-plot([x(z,t),[ksi(z,t),color=blue]],z=.83,t=0..0.3);

You will see an error message if you try t=0..0.4 or higher.
That probably is due to a singularity because of the denominator (1-Bi*(1-1/ksi(z,t))).
You can increase the upper t-limit slightly by using the optional argument timestep=0.001 in pdsolve.
So as I see it you have a problem if your interest is in t = 1..60.

This sure doesn't look much like Maple code to me. Where did you get it?

@Carl Love Thanks for the comment. I should have just said that the inspiration came from dsolve/numeric.
It is clearly convenient for the user of dsolve/numeric/listprocedure that it returns a list of equations whose left hand sides tell what the procedures on the right are computing.
Similarly dsolve/numeric/procedurelist returns a procedure which returns equations with the same left hand sides.
In general though this doesn't seem to be a good idea.

It would clearly be silly to use `convert/listprocedure` to the output from dsolve/numeric/procedurelist since dsolve/numeric/listprocedure is available and works fine.
However, for the fun of it we could try.
Example.
ivp:={diff(x(t),t,t)+x(t)=0,x(0)=0,D(x)(0)=1};
res:=dsolve(ivp,numeric);
#Output from numeric input to res is a list of equations
p:=proc(tt) subs(res(tt),[x(t),diff(x(t),t)]) end proc;
#Now p is a procedurelist as I'm using the term.
convert(p,listprocedure,2);
p1(Pi/2);
p2(Pi/2);
op(4,eval(p));
resL:=dsolve(ivp,numeric,output=listprocedure);
X0,X1:=op(subs(resL,[x(t),diff(x(t),t)]));
X0 and X1 should behave like p1 and p2.
X0(Pi/2);
X1(Pi/2);
####
#The type check numeric I have changed to complexcons since I posted the worksheet.

I have no problem. I just checked.
I'm using Windows 7 and Maple 17.02 Standard, worksheet mode, 1D input.