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.

@MarkusMPG On clicking on your link I get an error 404, no file etc.

@rwooduk You clearly cannot plot a graph of a function of 3 variables, for that you need 4 dimensions.
But you can fix one of the variables, e.g. k__z.
You need to remove &+-.

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);

A math program like Maple cannot help you set up the problem. What programs like Maple are designed for is solving the problem once it is set up.
After having set up the problem in mathematical terms you also have to feed some code to Maple. There you might get some help from this site if you encounter problems.