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

@rwooduk Could you upload a worksheet instead of giving a link to an image?
Use the thick green arrow in the icon panel of the editor.

@aylin 

I tried with the values of A[1] and A[2] you supplied. I took the endpoint to be 5 rather than 100, just to see if I could make it work.
I made the following changes:

Digits:=23:
#######

N:=5: #Endpoint, was 100
x0[1]:=5.5556*10^7; x0[2]:=1.1111*10^7;x0[3]:=6.3096*10^9;psiN[1]:=0;psiN[2]:=0;psiN[3]:=0;
#Using scaled variables y = x/x0 where x0 is the initial value:
PDEtools:-dchange({seq(x[i](t)=x0[i]*y[i](t),i=1..3)},{ode},[seq(y[i](t),i=1..3)]);
sys:=solve(%,{seq(diff(y[i](t),t),i=1..3),seq(diff(psi[i](t),t),i=1..3)});
icsy :={ y[1](0)=1, y[2](0)=1,y[3](0)=1,psi[1](N)=0,psi[2](N)=0,psi[3](N)=0};
res:=dsolve(sys union icsy,numeric,maxmesh=1024);

plots:-odeplot(res,[seq([t,y[i](t)],i=1..3)],0..5);
plots:-odeplot(res,[seq([t,psi[i](t)],i=1..3)],0..5);
plots:-odeplot(res,[seq([t,psi[i](t)],i=2..3)],0..5);

This worked, but I'm very sceptical about the result.
I tried N=10, but then dsolve went on for more than 10 minutes and then the error message:
Error, (in dsolve/numeric/BVPSolve) unable to achieve requested accuracy of 0.1e-5 with maximum 1024 point mesh (was able to get 54.), consider increasing `maxmesh` or using larger `abserr`

The worksheet you uploaded shows an error statement. None such occurs when I execute your worksheet unaltered.
This was discussed in
http://www.mapleprimes.com/questions/200393-How-To-Solve-4th-Order-Ode#comment201650

What are the values for A[1] and A[2] ?

@Jeroen In the statement
unapply~(list,x); #Notice the tilde ~

unapply is applied to each element of the list.

You could use map instead:

map(unapply,list,x);

Elementwise operation was introduced in Maple 13. See
?elementwise

@J4James Didn't you mean "I suspect that there is no exact solution"  (meaning no closed form) ?

@bails123 Did you use the thick green arrow in the icon panel? It's the last one in the second row.

If for some reason that proves impossible could you give us the Maple syntax from your worksheet as text?

There is no file attached.

@96mhaynes 
P must be negative since otherwise f''(x) > 0 for all x > 0 for which the solution exists, in which case f would never attain the value 0.
In fact we must demand at least initially that f''(0) < 0. Thus we at least have
eval(rhs(ode),{x=0,f(x)=1})<0;
which returns P < -1/20.

For guess work let us assume that f''(x) < 0 for 0 < x < a. Thus f'(x) > f'(a) = -1/sqrt(3) for 0 < x < a.
Integrating this inequality from x = 0 to x = a we get
f(a) - f(0) > (-1/sqrt(3))*a.
Thus by f(a) = 0 and f(0) = 1 we get a > sqrt(3).
This gives us some idea about the size of a; remember we don't know that f''(x) < 0 on the whole interval 0..a.

The guessing values for a and P in
sol:=fsolve([p1,p2],[1.9,-.6]);
conforms to the estimates for a and P given above.

Remark.
As it turns out our solution does indeed satisfy f''(x) < 0 for 0 < x < a:
plots:-odeplot(res,[x,eval(rhs(ode),P=sol[2])],0..10);

Additional remarks.
The maximal interval of definition of the solution found has upper bound something like 10.877911. Try
plots:-odeplot(res,[[x,f(x)],[x,diff(f(x),x)]],0..11);

If you are new to using the parameters option in dsolve/numeric you must realize that by doing
res(parameters=[-.5678]);
you are changing the parameter P in the ode and thus the solution. Plotting after a change in P gives a different graph (solution) than before. Notice that the procedure p sets the parameter and by using fsolve on p1,p2 many changes are made. The last value (when fsolve succeeds) corresponds to the one you want.