What algorithm is used in Matlab's pade function? I try to find out for quite some time now, how Matlab implements the calculation of Padé Approximants using its symbolic pade function. (the code of is buried in a compiled mex-file)
I compared its output with the "direct calculation" via coefficient comparison and solving the resulting linear systems of equations. I used floating point, variable precision arithmetic and pure symbolic variables and even in the last case Matlab's implementation achieves far better results for orders $> 7$.
So there must be a more sophisticated algorithm involved, but the documentation remains silent. The only tiny bit of information I found was the "Euclidian Algorithm" - but in my understanding that does not directly lead to the Padé-Approximant, just to the greatest common divisor of two numbers, am I right?
I contacted the Matlab Support and received the answer, that if no source is given the algorithm would be their own development. I actually can't believe that they developed an algorithm from scratch, if there are books full of alternatives (e.g. Baker, see below), or do they? At least their algorithm should be based on something?
Baker, George A.jun.; Graves-Morris, Peter, Padé approximants., Encyclopedia of Mathematics and Its Applications. 59. Cambridge: Cambridge Univ. Press. xiv, 746 p. (1996). ZBL0923.41001.
Can you give me some hints on how they could have implemented the algorithm? What approach for calculating Padé approximants could be based on the mentioned "Euclidian Algorithm"? How could I proceed from there?
I actually could just use their implementation, but I feel for any publication one should know how it is done.

Citing myself on https://en.wikipedia.org/wiki/Talk:Pad%C3%A9_approximant, one of the classical ways to compute Padé approximants uses the extended Euclidean algorithm. To compute the $[m/n]$ approximant, set 
\begin{align}
&&r_0&=x^{m+n+1},&v_0&=0,\\
&&r_1&=f(x),     &v_1&=1,\\
q_1&=r_0\text{ div }r_1,& r_2&=r_0-q_1r_1,& v_2&=v_0-q_1v_1\\
&\ \ \vdots\\
q_k&=r_{k-1}\text{ div }r_k,&r_{k+1}&=r_{k-1}-q_kr_k,&v_{k+1}&=v_{k-1}-q_kv_k,
\\
&\ \ \vdots
\end{align}
where $f(x)$ is the Taylor polynomial of degree $m+n$. Then in every step
$$
v_kf(x)=r_k\mod x^{m+n+1}
$$
which means that $\frac{r_k}{v_k}$ is the $[\deg r_k/\deg v_k]$ Padé approximant. For most $f$ the quotients $q_k$ will all have degree $1$ so that $\deg v_k=k-1$ and $\deg r_k=m+n+1-k$ so that in step $k=n+1$ one obtains the $[m/n]$ approximant.
