(Experimental) Can it be shown that this extension of the secant-interpolation has quadratic convergence?

Background: I needed some efficient but simple interpolation-methods aside of Newton's iteration because I want to have it in contexts, where the derivative of a function is not always known. So an obvious choice after bisection is regula-falsi and/or secant method. Now after a short course in wikipedia there it is mentioned, that the secant-method has convergence rate of about 1.6 (when it begins to converge), and a parabolic method named after Müller (using a quadratic interpolation-formula instead of the linear secant-formula) is a slight improvement.
Well, I thought that it would be interesting to extend this over the quadratic to cubic and finally dynamically increasing order, where the newly interpolated values and all old ones serve as a basis for the new estimate.

Rationale:
Assume some function $\small f(x)$; for a simple example I used $\small f(x) = exp(x)$. I want to find the coordinate $\small f(x_z)=z=5$ Meaningful bounds for the initial interval are for instance $\small X_2 = [1,3]$

I compute the vandermonde-matrix $\small Y_2$ as $\small Y_2=\underset{r=1..2, c=1..2}{\operatorname{matrix}} \left[ f(X_2[r])^c \right]$ ; then $\ [z,z^2] \cdot Y_2^{-1} \cdot X_2 = x_3$ gives an estimate for $\small x_z$ in the value of $\small x_3$ .

Now I append the new value $\small x_3$ at the vector $\small X_2$ to make it $\small X_3 = \operatorname{concat}(X_2,x_3)$. I create the quadratic vandermondematrix $\small Y_3$ in the obvious way and iterate with increased dimensions.(The secant-method analogue were simply not to increase the dimension but to put the new estimate as new first or second component into the X-vector, depending on the sign of $\small z - f(x_k)$ )

Empirically, if I'm near enough the true root, I get nearly exactly quadratic convergence, which is at least an improvement over the Müller's-method. Since the quadratic convergence with the Newton-method was in the wiki-article put as somehow a limit (for the cost of need of the derivative), this is then also interesting.
However, for not-so-well chosen initial points it seems at least as possibly erratic as again the Newton-method (as a simple example I just tried the function $\small f(x)=exp(x)-4*x^2$ having some quadratic disturbance over the exponential). And, since it requires inversion of increasing sizes of vandermonde-matrices with the additional problem, that with the increase of iterations the rows become massively linear depending, it is surely not a method for the practical use.

But for the principal understanding of that matter my questions are:
1) how could I attack the problem to show/argue/prove, that this process has in fact such a rate of convergence?
2) Is that quadratic convergence really a limit here? I hoped, that the dynamic of the method would as well allow a dynamic increasing of the rate of convergence.

Example: (iterations downwards; new estimates just appended to matrix, rightmost column log of absolute error)

$\qquad \small \begin{matrix} k & x_k=estim(x_z) & f(x_k) & \ln(|z-f(x_k)|) \\ 1 & 0 & 1.00000000000 & 1.38629436112 & \text{initial low x} \\ 2 & 5 & 148.413159103 & 4.96572968897 & \text{initial high x} \\ \hline \\ 3 & 0.135673098126 & 1.14530742925 & 1.34929125565 \\ 4 & 3.64058189042 & 38.1140084505 & 3.49995640888 \\ 5 & 3.30992901642 & 27.3831816466 & 3.10830985751 \\ 6 & 2.93800666159 & 18.8781781845 & 2.63031769174 \\ 7 & 2.53192098262 & 12.5776443810 & 2.02520238367 \\ 8 & 2.12577394716 & 8.37938017666 & 1.21769231299 \\ 9 & 1.79911612022 & 6.04430266580 & 0.0433493572149 \\ 10 & 1.63882927856 & 5.14913777620 & -1.90288472781 \\ 11 & 1.61020090108 & 5.00381639897 & -5.56844797984 & \text{convergence-rate quadratic}\\ 12 & 1.60943843944 & 5.00000263503 & -12.8466176558 \\ 13 & 1.60943791243 & 5.00000000000 & -27.3915350585 \\ 14 & 1.60943791243 & 5.00000000000 & -56.4731042545 \\ 15 & 1.60943791243 & 5.00000000000 & -114.629475172 \end{matrix}$

-
The standard reference on issues like these is Traub's book. Interpolatory methods like secant and Muller have a nonintegral convergence rate $p$, where $p$ is the positive root of the polynomial $p^{n+1}-p^n-\dots-1$ and $n$ is the order of the interpolating polynomial being used (e.g. $n=1$ for secant, $n=2$ for Muller). –  Ｊ. Ｍ. Oct 18 '11 at 16:08
Of course, the bigger the $n$, the nearer $p$ is to $2$... –  Ｊ. Ｍ. Oct 18 '11 at 16:12
"since it requires inversion of increasing sizes of vandermonde-matrices with the additional problem" - if you use Björck-Pereyra, the effort of solving Vandermonde equations is vastly reduced. –  Ｊ. Ｍ. Oct 18 '11 at 16:15
Ah, well, the book of Traub seems to be an exquisite resource for this; just skimmed a bit through two immediately relevant looking chapters. Thank you very much; I hope I'll find the book here. –  Gottfried Helms Oct 18 '11 at 16:40
Right; in that case, I insist that you look at Björck-Pereyra (which is actually divided differences in disguise); less effort is needed ($O(n^2)$ vs. LU's $O(n^3)$), and the solutions tend to be more accurate. –  Ｊ. Ｍ. Oct 18 '11 at 23:32