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

@Carl Love In 2D input space is interpreted as multiplication. One of the very good reasons to stay away from 2D input.

@96mhaynes fsolve is given two arguments: The first the list of 4 procedures each accepting 4 arguments a,b,f0,omega. The second argument to fsolve is a list of starting values for a,b,f0,omega. Good starting values may require some exploration and luck.

@96mhaynes You talked about a and b that is why I asked for two K's.

You must mean K = 1/sqrt(3) and not -1/sqrt(3).
Otherwise f'(a) would be positive since you said that f'(a) = -K.

I posted a solution with only the a-part.

@96mhaynes What are the values of the two K's?
In my simple example they were 1 and -1.

@96mhaynes If f has only 2 zeros a and b, then f´(a) + K =0 and f'(b) + K = 0 cannot both hold unless K = 0 (assuming smoothness of f).
Or do you use different values of K for a and b, one negative and the other positive?

@sarunas 
evalf(arctan(519/520));
convert(%,degrees);
evalf(%);
#Or in one step:
evalf(convert(arctan(519/520),degrees));