@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.

@rwooduk I don't see why it doesn't work.
Try downloading my worksheet:

MaplePrimes13-12-11p.mw

If this won't work either, you sure have a problem. Otherwise it could be a 2D-input problem. I stay away from 2D-input.

Indeed, it is the usual 2D-problem: In my code in the comment above I had a not very visible space between plot3d and the parenthesis. In 2D this is (stupidly) interpreted as a multiplication sign.
Try copying and pasting again exactly as you did. Then just below try
whattype(%);
You should find that the answer is `*`, which means that the expression last returned is a product.
Now try tho remove the space between plot3d and the parenthesis. Hit Enter. You should get a nice looking graph.

In 1D-input (aka Maple input) multiplication signs have to be inserted if you want to multiply.
Some people like 2D input and it is the default. However, you can easily change this behavior (and go back again):
Go to Tools/Options/Display/Input display and pick Maple Notation.
If you also want the worksheet interface (I do) then go on (still in Tools/Options) to
Interface/Default format for new worksheets and pick Worksheet.
Now click Apply Globally (no danger you can always go back).
Personally I don't like equation labels either. They can be unchecked under Display.

@samiyare This is an odd problem since it involves a system where a condition is given on the highest derivative appearing in one of the unknowns, in this case T.
Try
subs(eta=0,bcs1,convert(eq2=0,D));
and you will see that if phi[w] is real then T(0)=0.

And furthermore:
sol:=eval(dsolve({eq3=0,bcs1},phi(eta)),T(0)=0);
eq:=simplify(eval(eq2=0,sol));
res_T:=dsolve({eq,T(0)=0},T(eta));

You see that T(eta) = 0 for all eta so the problem has no solution satisfying also D(T)(0)=1.

@rwooduk Did you execute all the lines?

restart;
E:= (k__x, k__y, k__z)->
    gamma*sqrt(1+4*cos(3/2*k__y*a__cc)*cos(sqrt(3)/2*k__x*a__cc) +
               4*cos(sqrt(3)/2*k__x*a__cc)^2);

E:= subs([gamma= 3*1.6e-19, a__cc= 0.142e-9], eval(E));
a:=2e10:
plot3d ([E(k__x,k__y,a),-E(k__x,k__y,a)],k__x=-a .. a,k__y=-a .. a,shading=zhue);

@MarkusMPG It is hardly surprising that there are unevaluated integrals:

var:={p[c](1), p[c](2), p[c](3), p[c](4), p[c](5), p[c](6)};
remove(has,indets({odesys},function),var union {diff});
nops(%);
#12 unknown functions.

indets(Result,specfunc(anything,int));
nops(%);
#18 integrals are unevaluated for very good reasons: the integrands are not fully known.

@aylin I get the same error as before. Maybe it is a MaplePrimes editor problem.

@I_Love_LSK I would try another initialpoint.
However, I just checked in Maple 12. I get no such error when I do as described in the addendum to my last comment. The main points were:
1. scale
2. remove obvious (?) bad points.

xm:=max(op(map(abs,x11)));
ym:=max(op(map(abs,subsop(10=NULL,y11))));
zm:=max(op(map(abs,subsop(21=NULL,z11))));
#and then define ans:
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)/xm^2+add((Y(tim[i])-y11[i])^2,i in {$1..N} minus {10})/ym^2+add((Z(tim[i])-z11[i])^2,i in {$1..N} minus {21})/zm^2
 end proc;
Optimization:-Minimize(ans,initialpoint=[.001,.002,.003,.001,.002,.003,.001,.002,.003],method=nonlinearsimplex);
The sum of weighted squares reported is 202.4.
That can be improved by continuing with the result wec just got:
Optimization:-Minimize(ans,initialpoint=convert(op(2,%),list), method=nonlinearsimplex);
You can try the very same code line one more time. No improvement over the previous which had sum of squares 166.8.
Trying as I suggested as initialpoint the values I gave above (found in Maple 17) I get no improvement, but still a sum of squares = 6.7.