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

@aylin I get the following error when I click on either of the links:

404 - File or directory not found.
The resource you are looking for might have been removed, had its name changed, or is temporarily unavailable.

@aylin Is the paper available on the internet?

@Gaia Try using method=nonlinearsimplex:

Optimization:-Minimize(ans,initialpoint=[.001,.002,.003,.001,.002,.003,.001,.002,.003],method=nonlinearsimplex,evaluationlimit=1000);

That works in Maple 12. Also try with different initial points.
But as I said the result is not good.

You may try leaving out data point number 10 from y11 and number 21 from z11. This implies redefining 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)+add((Y(tim[i])-y11[i])^2,i in {$1..N} minus {10})+add((Z(tim[i])-z11[i])^2,i in {$1..N} minus {21})
 end proc;
Although that reduces the sum of squares by a factor 2*10^(-4) you still have a huge sum of squares.


Addendum. Because of the large differences between the x, y, and z values it seems reasonable to do some scaling. An easy way to do this is to apply weights to the sum of squares.
Leaving out data point 10 in y11 and 21 in z11 you can do:
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;

Then proceeding as before. After a while I got the somewhat decent (preliminary) result
[.196366324192122, 0.713565452459781e-2, 0.981673609018648e-1, -1.17885008868878, -.320685207261402, 1.50254628730238, -.596289885795190, 0.773099833270623e-2, -.323552878830983]
with the sum of weighted squares = 6.7.
This was done in Maple 17, but you may try in Maple 12 perhaps simply starting with the list values that I just gave.

@Carl Love Yes, I was too hasty!
I shall insert a note instead of deleting my response in view of the many times people have asked for all solutions.

Note added: For the univariate case

acer and Erik Postma has posted procedures a couple of years ago using NextZero:

http://www.mapleprimes.com/questions/123970-Finding-All-Roots-For-An-ODE-In-A-Given-Interval

http://www.mapleprimes.com/questions/123476-Seeking-A-Good-Praxis-In-Solving-Equations

They may well be worth looking at.

@Markiyan Hirnyk It looks as if the star is different from the centered dot, but that must be an artifact of the 2D input.
For the sake of clarity all 3 are the same:

(sqrt(M+m)*sqrt(M-m)*sqrt(M+m))*sqrt(M-m);