The substitution $z=y'$ is natural, so let us please then
forget about the posted $y$.
We will not use $z$. Let us use $y$ "back" instead of $z$ in the sequel.
So we solve (instead) around $x=0$:
$$ y'' ' = y\ y' \ /\ (1+x)\ .$$
Reusing the above code...
R10.<a0,a1,a2,a3,a4,a5,a6,a7,a8,a9,a10> = \
QQ['a0,a1,a2,a3,a4,a5,a6,a7,a8,a9,a10']
R.<x> = PowerSeriesRing(R10)
y = a0 + a1*x + a2*x^2 + a3*x^3 + a4*x^4 + a5*x^5 + \
a6*x^6 + a7*x^7 + a8*x^8 + a9*x^9 + a10*x^10 + O(x^11)
y1 = y .derivative()
y2 = y1.derivative()
y3 = y2.derivative()
f = (1+x)*y3 - y*y1
# J = ideal(f) # --> leads to an error
f.list()
... we get the coefficients of $f$:
[-a0*a1 + 6*a3,
-a1^2 - 2*a0*a2 + 6*a3 + 24*a4,
-3*a1*a2 - 3*a0*a3 + 24*a4 + 60*a5,
-2*a2^2 - 4*a1*a3 - 4*a0*a4 + 60*a5 + 120*a6,
-5*a2*a3 - 5*a1*a4 - 5*a0*a5 + 120*a6 + 210*a7,
-3*a3^2 - 6*a2*a4 - 6*a1*a5 - 6*a0*a6 + 210*a7 + 336*a8,
-7*a3*a4 - 7*a2*a5 - 7*a1*a6 - 7*a0*a7 + 336*a8 + 504*a9,
-4*a4^2 - 8*a3*a5 - 8*a2*a6 - 8*a1*a7 - 8*a0*a8 + 504*a9 + 720*a10]
We give or consider a0, a1, a2 as given.
Solving by the Cauchy-Kovalevskaya "procedure", we have to determine
recursively
a3,
a4,
a5,
a6,
a7,
a8,
a9,
a10,
and so on, starting from a0, a1, a2, by making each line evaluate to $0$.
We get manually
$$ a_3 = \frac 16 a_0a_1\ , $$
$$ a_4 = \frac 1{24}( a_1^2 + 2a_0a_2 + a_0a_1 )\ , $$
and so on.
I think it is not a good idea to isolate in this way all these equations, then start solving.
Instead, let us write the given DE in the form
$$ y'' ' = y\ y'\ (1-x+x^2-x^3+\dots)\ , $$
and consider it formally in the ring $\mathbb Q[[x]]$. Then the equations in $a_3,a_4,\dots$ come already
separated, and we can simply write the code to get the first some $100$ coefficients.
(I am not willing to write optimal code, time...)
To have a precise situation, let us better use values instead of the symbols $a_0,a_1, a_2$.
For instance $a_0=0$, $a_1=1$, $a_2=2$. Then...
a = [0,1,2] # use other initial conditions, next time
PSR.<x> = PowerSeriesRing( QQ, 100 )
C = 1 / PSR(1+x)
for k in [ 3..99 ]:
print k,
# trying to get a(k) from y = sum( [ a(k) * x^k for k in [0..ooh] ] )
yk = sum( [ a[j] * x^j for j in range(k) ] ) + O(x^k)
ak = 1/(k*(k-1)*(k-2)) * ( yk * diff( yk, x ) * ( C+O(x^k) ) ).list()[k-3]
print ak
a.append( ak )
y = PSR( sum( [ a[k] * x^k for k in [0..99] ] ) )
This computes some values.
For instance, the first some $30$ entries are:
sage: a[:30]
[0,
1,
2,
0,
1/24,
1/12,
1/40,
-67/5040,
13/1152,
-7/1440,
4621/1209600,
-38399/13305600,
11287/4838400,
-12737/6988800,
2351057/1614412800,
-171374993/145297152000,
1692465871/1743565824000,
-115101113/142502976000,
32202455561/47424990412800,
-4213094929/7313918976000,
22201339257697/45053740892160000,
-267905698350809/630752372490240000,
5115211318062343/13876552194785280000,
-8562348061276901/26596725040005120000,
928302494123858617/3282795776366346240000,
-7613558036179561/30490363247984640000,
1985785680757233642821/8962032469480125235200000,
-590031815878576041977/2987344156493375078400000,
33273819697787429364371/188202681859082629939200000,
-3468510498746975345435917/21831511095653585072947200000]
(These are exactly computed.)
Let us check the solution modulo $O(x^{100})$:
# : diff( y,x,3 ) - y * diff(y,x) / (1+x)
# the above is a mess with terms in degrees 97, 98, 99, but..
sage: diff( y,x,3 ) - y * diff(y,x) / (1+x) + O( x^97 )
O(x^97)
The radius of convergence is in this case probably one, one can try...
for k in [1..99]:
ak = a[k]
print ( "a[%2d] ~ %13.10f and a[%2d]^(1/%2d) ~ %8.6f"
% ( k, RR(ak), k, k, abs(RR(ak))^(1/k) ) )
and the first and last lines shown are:
a[ 1] ~ 1.0000000000 and a[ 1]^(1/ 1) ~ 1.000000
a[ 2] ~ 2.0000000000 and a[ 2]^(1/ 2) ~ 1.414214
a[ 3] ~ 0.0000000000 and a[ 3]^(1/ 3) ~ 0.000000
a[ 4] ~ 0.0416666667 and a[ 4]^(1/ 4) ~ 0.451801
a[ 5] ~ 0.0833333333 and a[ 5]^(1/ 5) ~ 0.608364
a[ 6] ~ 0.0250000000 and a[ 6]^(1/ 6) ~ 0.540742
a[ 7] ~ -0.0132936508 and a[ 7]^(1/ 7) ~ 0.539447
a[ 8] ~ 0.0112847222 and a[ 8]^(1/ 8) ~ 0.570902
a[ 9] ~ -0.0048611111 and a[ 9]^(1/ 9) ~ 0.553313
:::::::::::::::::::::::::::::::::::::::::::::::::
a[90] ~ 0.0000051223 and a[90]^(1/90) ~ 0.873406
a[91] ~ -0.0000049542 and a[91]^(1/91) ~ 0.874386
a[92] ~ 0.0000047934 and a[92]^(1/92) ~ 0.875348
a[93] ~ -0.0000046394 and a[93]^(1/93) ~ 0.876295
a[94] ~ 0.0000044920 and a[94]^(1/94) ~ 0.877225
a[95] ~ -0.0000043507 and a[95]^(1/95) ~ 0.878140
a[96] ~ 0.0000042154 and a[96]^(1/96) ~ 0.879040
a[97] ~ -0.0000040855 and a[97]^(1/97) ~ 0.879925
a[98] ~ 0.0000039610 and a[98]^(1/98) ~ 0.880796
a[99] ~ -0.0000038414 and a[99]^(1/99) ~ 0.881653
I have maybe to stop here with sage sample code.
About the convergence radius...
The recurrence relation is thus of the shape $a_k=$ a sum of
monomials of order two in the previous variables, with coefficients that can be
recursively "controlled". With some work, one may (??) maybe show that the radius of convergence is one.
(The one is my quick guess, based on the fact that the limit $r=\lim |a_k|^{1/k}$ may / should exist, and
then we expect $r\sim r^2$. Minoration / majoration for the "polynomial part in $k$" is constricted by the power $1/k$. But it is just a guess.)
Note: The use of solve_laplace really wants to perform a Laplace transform. But which is the first step
when we try to apply it on the right side of the given DE ?!
Note: The answer is not given mot-a-mot to the posted (relatively general) question, but in the sense of the progress. (Which is subjective.) I think, sample code serves good as an answer.
