prove conjecture; the limit of iterating is $\sqrt{z^2 - 2}$ $$\lim_{n \to \infty} f_n(x)=x-\frac{1}{nx}\;\;\; g(x)=f_n^{on}(x)$$
The conjecture is for values of $|x|>\sqrt{2}$: $g(x) = \sqrt{z^2 - 2}$
This question comes from another matstack question/answer.
If $n=2^m$, then convergence is much quicker by starting with $f_n(x)=x-\frac{1}{nx}$ and then iterating the Taylor/Laurent series for $f(x) \mapsto f(f(x))$ m times.  I start by generating the formal Taylor series after moving the fixed point from infinity to zero, $1/f(\frac{1}{x})\;$.  Then we start iterating with $f_n$ 
$$f_n(x) = \frac{x}{1 - x^2/n} =  x + \frac{x^3}{n} + \frac{x^5}{n^2} + \frac{x^7}{n^3} + \frac{x^9}{n^4} + \frac{x^{11}}{n^5} ...$$
Using this speedup with $m=\log_2(n)$, one can iterate the Taylor series for $f(x) \mapsto f(f(x))$ m times, rather then iterating $f^{on}$, but the two are of course identical.  Its just that otherwise convergence is pretty slow, with n iterations to get accuracy to 1/n.
The formal Taylor series coefficients of $1/g(\frac{1}{x}) = \sqrt{z^2/(1-2z^2)}$ are:
$x + x^3 +
\frac{3        x^5    }{2     } +
\frac{5        x^7    }{2     } +
\frac{35       x^9    }{8     } +
\frac{63       x^{11} }{8     } +
\frac{231      x^{13} }{16    } +
\frac{429      x^{15} }{16    } +
\frac{6435     x^{17} }{128   } +
\frac{12155    x^{19} }{128   } +
\frac{46189    x^{21} }{256   } ...$
Empirically, for the $2^n$th iteration starting with $f_n$ in the limit above, the Taylor coefficients are accurate to approximately $O 2^{-n}$, so there is pretty good numerical computation evidence for the conjecture, but I have no idea how to prove it.
EDIT
I (https://math.stackexchange.com/users/39261/mick) will place a bounty for the following problem :
Basicly the inverse :
Suppose we are given $g(x)=\sqrt{x^2-2}$ and we are asked to find $f_n$ such that :
$$\lim_{n \to \infty} \;\;\; g(x)=f_n^{on}(x)$$
How do we solve such problems ??
EDIT
EDIT 2
As Sheldon's comment says $f_n(x)=\sqrt{x^2-\frac{2}{n}}$ is also a solution but I want to find the $f_n$ from the OP : $x - \frac{1}{nx}$ or another $f_n$ that is real-meromorphic on the entire complex plane. 
EDIT 2
 A: I'll solve a slightly more general problem with this post - that is, if $f_n(x)=x-\frac{1}{nx}$, what is the limit
$$g(x)=\lim_{n\rightarrow\infty}f^{n}_n(x).$$
To do this, we first change the form of $f_n$ to a conjugation of other functions. If we define:
$$f(x)=x-\frac{1}x$$
$$\alpha_n(x)=\sqrt{n}\cdot x$$
then
$$f_n=\alpha_n^{-1}\circ f \circ \alpha_n$$
The advantage of this is that we can rewrite, where $g_n=f_n^n$ that
$$g_n = \alpha_n^{-1}\circ f^n \circ \alpha$$
which is nice because it means we only have to study the iterates of $f$ to figure this out. However, this isn't too bad, since we only need tight bounds on $f$. In particular, consider the following equation:
$$\Delta f^n=\frac{-1}{f^n}$$
where we take $\Delta_n f_n$ to be defined as $f^{n+1}-f^n$. I write it this way, because next, we can make a leap to a related differential equation; in particular, let's a function $I_x(n)$ which will approximate $f_n$. In particular, we want the following equations to hold:
$$I_x(0)=x$$
$$I_x'(n)=\frac{-1}{I_x(n)}$$
This is a nice separable differential equation and has the solution (for positive $x$):
$$I_x(n)=\sqrt{x^2-2n}.$$
The precise connection here is that $f^n(x)$ is basically the result of applying Euler's (forward) method to the same differential equation, with step size $1$. However, notice that we have the identity:
$$\frac{I_{\sqrt{k}x}(kn)}{\sqrt{k}}=I_x(n)$$
which tells us that our differential equation scales in a certain way (see remarks at end if you want another way to see this) - but what's remarkable, is this means that $f^{kn}(\sqrt{k}x)$ uses Euler's method to approximate $I_{\sqrt{k}x}(n)$ with step size $1$ - but when we scale this result back by the above equation, we get that $\frac{f^{kn}(\sqrt{k}x)}{\sqrt{k}}$ is actually Euler's method, with step size $\frac{1}k$, approximating $I_x(n)$. As the differential equation is nicely behaved (i.e. Lipschitz continuous), we reach the result that, as we decrease the step size, we get
$$\lim_{k\rightarrow\infty}\frac{f^{kn}(\sqrt{k}x)}{\sqrt{k}}=I_{x}(n)$$
setting $n=1$ yields:
$$\lim_{k\rightarrow\infty}\frac{f^{k}(\sqrt{k}x)}{\sqrt{k}}=I_{x}(1)$$
and substituting
$$g(x)=\sqrt{x^2-2}$$
as desired.

To really be clear, it's wise to note that when we're talking about Euler's method, what we're really doing is taking the points given from the scaling operation $\frac{1}{\sqrt{k}}f^{kn}(\sqrt{k}x)$ - which are points $kn$ which are $(n,f)$ pairs:
$$\left(0,x\right),\left(\frac{1}{k},\frac{1}{\sqrt{k}}f(\sqrt{k}x)\right),\left(\frac{2}k,\frac{1}{\sqrt{k}}f^2(\sqrt{k}x)\right),\ldots$$
and saying, "Gee, that looks like Euler's method". Though, given the identity on $I_x(n)$ it is clear that this must work, it's good to try an example to be sure. Notice that the second point above ought to be equal to $x+\frac{1}kI_x'(0)$, as would be Euler's method. This ends up expanding as:
$$x+\frac{1}kI_x'(0)=\frac{1}{\sqrt{k}}f(\sqrt{k}x)$$
$$x-\frac{1}{kx}=\frac{1}{\sqrt{k}}(\sqrt{k}x-\frac{1}{\sqrt{k}x})$$
$$x-\frac{1}{kx}=1-\frac{1}{kx}$$
Since all pairs of points have this same basic form, we can prove algebraically (and inductively) that $\frac{f^{kn}(\sqrt{k}x)}{\sqrt{n}}$ is actually Euler's method approximating $I_x(n)$ with step size $\frac{1}k$ which does suffice to prove the claim. This warped my mind for an hour, but I think I'm okay now.

A particular generalization of this that we can make is that if we let
$$f_n(x)=x+\frac{\alpha}{2nx}$$
and, correspondingly
$$f(x)=x+\frac{\alpha}{2x}$$
then we have
$$\Delta f^n=\frac{-\alpha}{f^n}$$
and correspondingly
$$I_x'(n)=\frac{\alpha}{2I_x(n)}$$
which has the solution
$$I_x(n)=\sqrt{x^2+\alpha n}$$
and by the same argument, we prove that $$\lim_{n\rightarrow\infty}f_n^n(x)=\sqrt{x^2+\alpha}$$
A: Pick an $x > \sqrt{2}$ and $n > \dfrac{1}{x^2-2}$ and define a sequence by $x_0 = x$, $x_{k+1} = f_n(x_{k})$ for $k \ge 0$. 
Rearrange $x_{k+1} = f_n(x_k) = x_k - \dfrac{1}{nx_k}$, to get $x_k(x_k-x_{k+1}) = \dfrac{1}{n}$. 
Clearly, $x_0 \ge \sqrt{x^2}$. Now, suppose that $x_{m-1} \ge \sqrt{x^2-\dfrac{2(m-1)}{n}}$ for some $1 \le m \le n$. 
Then, $x_{m-1} \ge \sqrt{x^2-2} > \dfrac{1}{\sqrt{n}}$, and so, $x_m = x_{m-1} - \dfrac{1}{nx_{m-1}} > 0$. Furthermore, we have: 
$x^2-x_m^2 = \displaystyle\sum_{k = 0}^{m-1}(x_k^2-x_{k+1}^2) = \sum_{k = 0}^{m-1}\left[(2x_k-(x_k-x_{k+1}))(x_k-x_{k+1})\right]$ $= \displaystyle\sum_{k = 0}^{m-1}\left[2x_k(x_k-x_{k+1}) - (x_k-x_{k+1})^2\right] = \sum_{k = 0}^{m-1}\left[\dfrac{2}{n} - \dfrac{1}{n^2x_k^2}\right] = \dfrac{2m}{n} - \dfrac{1}{n^2}\sum_{k = 0}^{m-1}\dfrac{1}{x_k^2} \le \dfrac{2m}{n}$. 
Hence, $x_m^2 \ge x^2 - \dfrac{2m}{n}$. Then, since $x_m > 0$, we have $x_m \ge \sqrt{x^2-\dfrac{2m}{n}}$. 
So by induction, we have $x_n \ge \sqrt{x^2-2}$. (Note: all of that was necessary to ensure that $x_n > 0$ , otherwise, $x^2-x_n^2 \le 2$ doesn't imply that $x_n \ge \sqrt{x^2-2}$.)
Then, since $x^2-x_n^2 = 2 - \dfrac{1}{n^2}\displaystyle\sum_{k = 0}^{n-1}\dfrac{1}{x_k^2} \ge 2 - \dfrac{1}{nx_n^2} \ge 2 - \dfrac{1}{n(x^2-2)}$, we have $x_n^2 \le x^2-2 + \dfrac{1}{n(x^2-2)}$, i.e. $x_n \le \sqrt{x^2-2 + \dfrac{1}{n(x^2-2)}}$.
Therefore, $\sqrt{x^2-2} \le f_n^n(x) \le \sqrt{x^2-2 + \dfrac{1}{n(x^2-2)}}$ for all $n > \dfrac{1}{x^2-2}$. 
Taking the limit as $n \to \infty$ and squeezing yields $g(x) = \displaystyle\lim_{n \to \infty}f_n^n(x) = \sqrt{x^2-2}$, as desired. 
Note that the proof that $g(x) = \sqrt{x^2-2}$ for $x < -\sqrt{2}$ is similar.
A: @Mick comments:
As Sheldon's comment says $f_n(x)=\sqrt{x^2-\frac{2}{n}}$ is also a solution but I want to find the $f_n$ from the OP : $x - \frac{1}{nx}$ ...
The answer to Mick's comments might
be as follows, which would be another solution to the Op's question.  Here is a formal equation for the Abel function for $g(z)$
$$f(z)=z-\frac{1}{zn}$$
Let $\alpha_f(z)$ be the formal Abel function solution such that
$$\alpha_f(f(z)) = \alpha_f(z)+1\;\;\;g(z) = f(z)^{on}$$
Therefore $$
\alpha_g(z) = \frac{\alpha_f(z)}{n};\;\;\; \alpha_g(f(z))=\alpha_g(z)+\frac{1}{n};\;\;\; \alpha_g(g(z))=\alpha_g(z)+1$$
To calculate $\alpha_f(z)$, we use the method of Ecalle, which works
for parabolic fixed points.  First, we move the fixed point of $f_n$ from
infinity to zero.
   $$ f_n(z) = z-\frac{1}{nz} \mapsto \frac{1}{f_n(1/z)} $$
   $$ \frac{1}{f_n(1/z)} = \frac{z}{1 - z^2/n} =  z + \frac{z^3}{n} + \frac{z^5}{n^2} + \frac{z^7}{n^3} + \frac{z^9}{n^4} + \frac{z^{11}}{n^5} ...$$
Then using the attatched pari-gp program, one gets the assymptotic equation for the Abel function, $\alpha_f(z)$, and as noted, $\alpha_g(z)=\alpha_f(z)/n$.  The first few terms of the formal $\alpha_g(z)$ series are printed below. To get arbitrarily accurate convergence, one would iterate $f^{-1}(z)$ a few times, $\alpha_g(z)=\alpha_g(f^{o-m}(z))+\frac{m}{n}$
$$\alpha_f(z) = 
\frac{-n}{2z^2} +
\frac{\ln(z)}{2} +
\frac{-z^2}{8n} +
\frac{5z^4}{96n^2} +
\frac{-7z^6}{288n^3} +
\frac{-1z^8}{1280n^4} +
\frac{671z^{10}}{28800n^5} +
\frac{-9607z^{12}}{483840n^6} ... $$
$$\alpha_g(z) = 
\frac{-1}{2z^2} +
\frac{\ln(z)}{2n} +
\frac{-z^2}{8n^2} +
\frac{5z^4}{96n^3} +
\frac{-7z^6}{288n^4} +
\frac{-1z^8}{1280n^5} +
\frac{671z^{10}}{28800n^6} +
\frac{-9607z^{12}}{483840n^7} ... $$
$$ \alpha_g(g(z)) = \alpha_g(z)+1$$
For arbitrary fractional iterations of $g(z)$, for finite values of n, we make due with numerical methods.  Here, since we moved the fixed point to zero, $g(z) \mapsto \frac{1}{g(1/z)}$ 
$$g^{oz} = \alpha^{-1}(z)$$ 
$$g(z)=\alpha_g^{-1}(\alpha_g(z)+1)$$
Conveniently, as n gets arbitrarily large, the Abel function converges, which leads directly to the solution.  
$$\alpha(z)=\frac{-1}{2z^2};\;\;\; \alpha^{-1}(z) = -\sqrt{\frac{-1}{2z}}$$
And then with a little algebra
$$g(z)=\alpha_g^{-1}(\alpha_g(z)+1);\;=\sqrt{\frac{z^2}{1-2z^2}}$$
Finally, moving the fixed point back to infinity yields
$$g(z) \mapsto \frac{1}{g(1/z)} = \sqrt{z^2-2}$$
Here is the pari-gp program I used to calculate the formal Abelseries of $\alpha_g(z)$, which is based on Ecalle's method.  See Will Jagy's math overflow post
\ps 24
/* kn is just n, but I use n already so ... */
kx=0; forstep (n=0,12,1,kx=kx+x^(2*n+1)/kn^n);
abelseries(fz,n) = {
  local(i,z,ns,rem,logfzx);
  kabel=0;
  klog=0;
  kfz=fz; /* save fz */
  logfzx=Ser(log(fz/x));
  negterms=1;
  while (polcoeff(fz,negterms+1)==0,negterms++);
  print("terms with negative coeffients= "negterms);
  for (i=-negterms,n,
    if (i==0, klog=acoeff, kabel=kabel+acoeff*x^i);
    rem = Ser(subst(kabel,x,fz) - kabel + klog*logfzx - 1);
    z=polcoeff(rem,i+negterms);
    z=subst(z,acoeff,x);
    ns=-polcoeff(z,0)/polcoeff(z,1);
    kabel=subst(kabel,acoeff,ns);
    klog=subst(klog,acoeff,ns);
  );
  return([kabel,klog]);
}
kz = abelseries(kx,20);
kg = kz/kn;

