Often in physics problems, we derive some kind of function f(x) for some physical quantity and we try to find the approximate value of the function as some variable gets too big, or close to zero.
To do this formally, I take the limit of the function as normal. But what is the formal way to approximate the function at x much bigger than another variable? I mostly see it stated as $x>>y$. A simple example would be $|E(z)|=\dfrac{1}{4\pi\epsilon_0} 2\pi\sigma (1 - \dfrac{z}{\sqrt{z^2+R^2}})$ for the electrical field of a uniformly charged circular disk, with radius $R$, at distance $z$ from the center of the disk.
For $z \to 0$, the second term is equal to $0$.
For $z \to \infty$, $\lim_{z \to \infty}\dfrac{z}{\sqrt{z^2+R^2}} = \lim_{z \to \infty}\dfrac{z}{z\sqrt{1+(\dfrac{R}{Z})^2}} = \lim_{z \to \infty}\dfrac{1}{\sqrt{1+(\dfrac{R}{z})^2}} = 1$
But how do I go about formally calculating the second term for $z >> R$ and $z << R$? Do I have to express the function as $E(u)=E(z-R)$ and then take the limit as $u \to \infty$? In most cases, building up the $z >> R$ expression works fine without much trouble but is there a specific methodology I could follow?
My textbook provides the solution as follows:

which feels kind of random to me.
The problem here is that there isn't just one approximation that will work. Even though one likes to say things like "when $x$ is large, then $f(x) \approx g(x)$ (for some 'simpler' $g$)", that does not imply that the $g$ is unique. After all, it still is an approximation even "when $x$ is large", and that means we can think of better ones, which will necessarily have a different form.
I would venture to say this is more "art" than analytical. That said, we can try to put some theory to it. Consider the $x \rightarrow \infty$ case for simplicity. Let $f_j$ be an infinite sequence of functions such that we can form a very general infinite series
$$f(x) = f_1(x) + f_2(x) + f_3(x) + \cdots$$
where this equality is exact, for at least $x$ greater than some minimum value $x_0$. Suppose each function "grows slower" than the last in that $f_{j+1}(x) \in o(f_j(x))$. Then if we truncate the series after some number of terms, the error, represented by the terms remaining that we removed, will die off as $x$ grows.
However, there is virtually unlimited freedom in the form the $f_j$ may take. For example, we could just take $f_1 = f$, then for all other $j$, $f_j = x \mapsto 0$. In that case, the above property will hold, except it is not particularly useful: we can "approximate" $f$ to any order we like by $f$ itself. And that's where the "art" comes in: it lies is in choosing the most convenient sequence of "terms" $f_j$ to express $f$ as so we can then truncate it.
And one common, but not always applicable, way to do that is to use a Taylor or Laurent series expansion to generate a suitable sequence of $f_j$. In your problem, we express
$$z \mapsto \frac{1}{\sqrt{z^2 + R^2}}$$
by noting
$$\frac{1}{\sqrt{z^2 + R^2}} = \frac{1}{z} \left(1 + \frac{R^2}{z^2}\right)^{-1/2}$$
and then we can use the binomial series
$$(1 + z)^r = \sum_{n=0}^{\infty} \binom{r}{n} z^n$$
to get a Laurent series
$$\left(1 + \frac{R^2}{z^2}\right)^{-1/2} = \sum_{n=0}^{\infty} \binom{-1/2}{n} \frac{R^{2n}}{z^{2n}}$$
hence, multiplying by $\frac{1}{z}$,
$$\frac{1}{\sqrt{z^2 + R^2}} = \sum_{n=0}^{\infty} \binom{-1/2}{n} \frac{R^{2n}}{z^{2n + 1}} = \frac{1}{z} - \frac{1}{2} \frac{R^2}{z^3} + \frac{3}{8} \frac{R^4}{z^5} - \frac{5}{16} \frac{R^6}{z^7} + \cdots$$
we have
$$f_1(z) = \frac{1}{z}$$ $$f_2(z) = -\frac{1}{2} \frac{R^2}{z^3}$$ $$f_3(z) = \frac{3}{8} \frac{R^4}{z^5}$$ $$f_4(z) = -\frac{5}{16} \frac{R^6}{z^7}$$ $$\cdots$$
and you can show that, using the definition of $o$, $\frac{1}{z^3} \in o\left(\frac{1}{z}\right)$, $\frac{1}{z^5} \in o\left(\frac{1}{z^3}\right)$, etc. Hence we have a sequence that meets the theory's prescriptions. Truncating to 2 terms then gives you the result in the textbook.