@acer 

I did not do it by hand, but Maxima agrees.

Edit. Maybe  99.9999 because Derive agrees too (yes, I still use it sometimes!)

This reminds me of:

@williamdiao 

The only problem is that Quo is called with 2 or 3 arguments. So, e.g. use:

C[Quo] := proc(a,b,r) if nargs=3 then quo(a,b,z,r) else quo(a,b,z) fi end:

@Math111 

Then what is (should be) F''(0) ?
[the fact that a method provides a solution does not imply that it actually exists].

Edit.

RK-4 for a BC problem? Probably some shooting method.

The ODE has the following explicit form:

The denominator has big chances to be 0 for r in [0,1].

Note that:
|F'(r)| must be < 0.25 ==>  |F(r)| < 1 + 0.25 r

The usual way of presenting this would be:

Write the projection as a procedure. E.g. the projection onto the unit sphere is
P:=(x,y,z) -> (x/sqrt(x^2+y^2+z^2), y/sqrt(x^2+y^2+z^2), z/sqrt(x^2+y^2+z^2));

Define the curve. E.g.
C:=t -> (cos(t),sin(t),t);

Plot the projection:
plots:-spacecurve( [P(C(t))],t=0..8*Pi);
or
plots:-animate(plots:-spacecurve, [[P(C(t))],t=0..a],a=0..8*Pi );

Then your posts will have a larger audience.

@Simwar 

For rather small intervals you can check all the possible cases.

Once you have a plausible formula, you could try to compute the sum using Maple (with assume facility, maybe fixing some parameters). Maple is "competent" in computing this kind of sums.

But if you have derived correctly the formula, a comprehensive random test should be enough for your peace of mind.

@Ronan 

Yes, it's odd that the rotation matrix for rotate(Plot, u,v,w)  is Rz(-w).Ry(-v).Rx(-u).

@Ronan 

In this new version you should use:

Bdyprt1 := unapply('display'('point'([d[1][1], d[1][2], d[1][3]]), colour = red, symbolsize = 50, symbol = solidsphere), t):
Bdyprt1a := 'animate'('Bdyprt1', [t], t = 0 .. 8*Pi, frames = 100):

etc
(but your approach does not seem very "natural").

Unfortunately the situation seems to be worse in Windows for 3d vector graphics:

http://www.mapleprimes.com/questions/220217-Bad-Quality-When-Exporting-To-PDF#answer234206

@acer 

I agree!  But you must admit that it is not normal that a syntax accepted by  line is rejected by its twin brother points.

@acer 

The unpleasant fact is that absolutely similar procedures such as line and points have distinct evaluation rules which must be guessed by trial and error.

@roman_pearce 

Yes, I know the theorem. It applies usually to extend (y0,z0) to (x0,y0,z0) and here the extension is not always possible, even in very simple examples. (I have used the notations from the question).

But I have not seen an impossible extension from (z0) to (x0,y0,z0). Have you?

Thank you for the answer.

@Ronan 

point does some strange manipulations. Workaround:

ff:=unapply(display(POINTS([d[1][1], d[1][2], d[1][3]]), colour = red, symbolsize = 50, symbol = solidsphere), t);

@mehdi jafari 

The syntax problem is that (applied to your example):

fsolve({seq(eq||i,i=1..3)},{a1=-5..0,a2,a3});  #error
fsolve({seq(eq||i,i=1..3)},{a1=-5..0,a2=-infinity..infinity,a3=-infinity..infinity}); #ok

(Or set the ranges separately, after vars, as you did).

@williamdiao 

Have you read the answer+comments+example below?