Thank you very much jakubi, this is outstanding. Very useful to understand this Panayotounakos formula. I am also learning a great deal about how to use Maple's capabilities.

The discrepancies between the dsolve,numeric output and the Panayotounakos formula can get quite large, what does it mean?

Take the simple case f(x)=0. In this case, the ODE can be solved in closed form. You'd expect numerical inaccuracies to be minimal in a simple case such as this. Yet, the discrepancies do look rather large.

Consider the case F(t)=0 . I input the Panayotounakos formula a little differently from jakubi, but it amounts to the same: I define a function H(t), in terms of cos(t), sin(t), and Ci(t), and express the Panayotounakos formula in terms of H(t), and F(t). (The formulae could be made even cleaner by defining a function exp(-t)*F(t) and factorizing exp(-2*t), but I haven't bothered) Then I input the trigonometric formula for the cubic roots directly, rather than using implicitplot -- the graph that obtains in this way, for the case F(t)=exp(-t), is the same as the one jakubi has presented. I then looked at different rhs functions, and compared the dsolve,numeric output with the formula -- discrepancies seem to appear even far from the discontinuities.

Starting with the transformed Abel ODE in r(t):

de := (r(t)+1/3)*diff(r(t),t)+r(t)^2-4/3*r(t)-5/9 = 4*F(t)/exp(t); #Transformed Abel ODE

H := t-> ( ( t*sin(t) + cos(t) ) * Ci(t) + (cos(t))^(2) ) * ( 4*t*Ci(t) + cos(t) ) / ( (t*Ci(t))^(3) ); #Panayotounakos's function

P := (t,H,F) -> -7/9 + 1/12 * H * exp(-2*t) - 4/3 * F * exp(-t) ;

Q := (t,H,F) -> -10/27 + 1/24 * H * exp(-2*t) - 8/3 * F * exp(-t) ;

R := (P,Q,k)-> 2*(-P)^(1/2) * sin(( arcsin(Q/(P)^(2)*(-P)^(1/2)) +2*k*Pi)/3); #trigonometric formula for roots

K := 1: # Selecting one of the real branches of the cubic, K=0,1,2

rsol := (t,F)-> R(P(t,H(t),F),Q(t,H(t),F),K): rsol(t,F): #Panayotounakos's formula for the solution r(t)

The above completes the coding of the Panayotounakos formula for the solution (to the best of my understanding). I now take a function on the rhs of the ODE and compare the formula with dsolve,numeric, following jakubi's suggestions.

FF := t-> 0: #Test function for F(t)

rsol1 := t-> rsol(t,FF(t)): rsol1(t):

r1 := evalf(Re(rsol1(1.)));

r1 := 1.16059739086600

de1 := subs( F(t) = FF(t), de );

dsolve(de1); #In some cases the solution can be found explicitly

r(t) = -1/3, r(t) = 5/3 + exp(-t) _C1

de1sol := dsolve( { de1, r(1.) = r1 }, numeric );

p1 := plot( rsol1(t) , t=0.2..4*Pi, resolution=1000, color=green ):

p2 := odeplot( de1sol, 0.2..4*Pi, resolution=1000, color=blue ):

display({p1,p2}, view=[1..4*Pi,0..3]);

Anyhow, I need to think about all this some more. I managed to post a picture, so that's not bad going for one day.

Thanks jakubi for your great help!

Patrick.