Regarding Padé Approximation of Neumann Series I calculated the Padé Approximation of Neumann Series by hand, and then by Mathematica for different orders (from {0, 0} to some higher numbers), using the code below, in general:
neu = 1 + x + x^2 + x^3 + x^4 + x^5 + x^6 + x^7 + x^8;
PadeApproximant[neu, {x, 0, {1, 1}}]

I got exactly the same results. It must be good news, but according to some articles, I should get different results for different orders of this approximation. For any order higher than {1, 1}, such as {2, 1}, {2, 2}, {3, 2}, ..., the result is $\frac{1}{1-x}$ (which is of course the result for the summation of an infinite number of terms of Neumann series). I just need to be sure that these results are correct and it makes sense that after {1, 1}, no progress happens in this approximation at all. Thank you very much in advance.
 A: The Pade approximation is a rational function that best matches the polynomial. Since 1/(1-x) already matches the terms you provide, there's nothing additional for higher-order Pade to approximate.
Another way to think of it, is the coefficients of those higher order terms happen to be zero for this particular polynomial.
A: Since this question is posed in Mathematica format, I first answer it in Mathematica format.  Consider
pa = FullSimplify@PadeApproximant[f[x], {x, 0, {8, 1}}]
(* f[0] + (x*(56*(720*Derivative[1][f][0] + 
   6*x*(60*Derivative[2][f][0] + 20*x*Derivative[3][f][0] + 
   5*x^2*Derivative[4][f][0] + x^3*Derivative[5][f][0]) + 
   x^5*Derivative[6][f][0]) + 8*x^6*Derivative[7][f][0] - 
   (9*x^7*Derivative[8][f][0]^2)
   /(-9*Derivative[8][f][0] + x*Derivative[9][f][0])))/40320 *)

Now, for f[x] given by neu, the first eight Derivatives evaluated at x = 0 are equal to n!, and for higher derivatives by 0.  So, clearly the expression above cannot have a denominator proportional to (1 - x).  In fact, it is given by
Simplify[pa /. {f[0] -> 1, Derivative[n_][f][0] -> If[n <= 8, n!, 0]}]
(* 1 + x + x^2 + x^3 + x^4 + x^5 + x^6 + x^7 + x^8 *)

On the other hand, if neu were a ninth order series, then the result would be
Simplify[pa /. {f[0] -> 1, Derivative[n_][f][0] -> If[n <= 9, n!, 0]}]
(* 1/(1 - x) *)

More generally, 1/(1 - x) is obtained for other Pade Approximations of neu that do not involve 9th or higher order derivatives (and, of course, have at least first order denominators).  The results were obtained with Mathematica
$Version
(* "13.0.0 for Microsoft Windows (64-bit) (December 3, 2021)" *)

A verbal explanation of these results is that neu is equal to 1/(1 - x) near x = 0 to eight order, but not to higher order.  So, low order Pade approximations equal 1/(1 - x), but high order Pade approximations are better approximated by neu itself.  I hope this is helpful.
