This is more a comment, addressing rather the optimal setting to get many iterations in time. (It became an answer since there is no fit in space.)
The short answert to the question is to take j( x, 0 ), which is still has the parent Multivariate Power Series Ring in x, y over Rational Field. See the detailed answer of tmonteil .
The long answer following is hopefully not so long.
The limit is of course known, we solve in the fix point equation $y=f(x,y)$ for $y$. I will take for $f$ the version with the plus sign, $f(x,y)=x^3+xy^2$. The solution to this equation, rewritten as $xy^2 -y +x^3=0$ is formally
$$ y =frac 1{2x}( 1\pm \sqrt{ 1-4x^4} )\ ,$$
which can also be easily implemented. The nonsingular solution is obtained for the choice of the minus sign above. The following is the code for the iterations
$$y_0=y\ , \qquad y_{n+1}=f(x,y_n)\ ,$$
showing from them the terms that are established from the expected formal limit, plus a four in the valuation. (Next, unestablished term is also shown.) It is important in code to have a truncation when defining f, else the 8.th iteration is still hard to compute in a glance.
The following code, written in the spirit of the posted question shows the linear convergence:
Code:
PREC = 50
STEPS = 7
R.<x, y> = PowerSeriesRing( QQ, default_prec=PREC )
f = x^3 + x*y^2 + R.O( PREC )
LIM = ( 1 - exp( log( 1-4*x^4 ) / 2 ) ) / 2 / x + R.O( PREC )
print "Expected limit:\nLIM = %s" % LIM
fn = y # f0 - the start of the iteration process
for n in [ 1..STEPS ]:
fn = f( x, fn )
yn = fn( x, 0 )
val = ( LIM - yn ).valuation()
print "f%s = %s" % ( n, yn + R.O( val+4 ) )
Results:
Expected limit:
LIM = x^3 + x^7 + 2*x^11 + 5*x^15 + 14*x^19 + 42*x^23 + 132*x^27 + 429*x^31 + 1430*x^35 + 4862*x^39 + 16796*x^43 + 58786*x^47 + O(x, y)^49
f1 = x^3 + O(x, y)^11
f2 = x^3 + x^7 + O(x, y)^15
f3 = x^3 + x^7 + 2*x^11 + x^15 + O(x, y)^19
f4 = x^3 + x^7 + 2*x^11 + 5*x^15 + 6*x^19 + O(x, y)^23
f5 = x^3 + x^7 + 2*x^11 + 5*x^15 + 14*x^19 + 26*x^23 + O(x, y)^27
f6 = x^3 + x^7 + 2*x^11 + 5*x^15 + 14*x^19 + 42*x^23 + 100*x^27 + O(x, y)^31
f7 = x^3 + x^7 + 2*x^11 + 5*x^15 + 14*x^19 + 42*x^23 + 132*x^27 + 365*x^31 + O(x, y)^35
Only seven iterations that fit in the line were done and shown here.
