You have several choices about how to handle the limits of integration.
The most straightforward thing to do is to observe that for any given $x$, the $y$-coordinate of a point in the disk must lie between $-\sqrt{1-x^2}$ and $+\sqrt{1-x^2}$, since $y$ outside this range makes $x^2+y^2 > 1$. So the integral is:
$$\int_{x=-1}^{x=1}\int_{y={-\sqrt{1-x^2}}}^{y=\sqrt{1-x^2}} \sigma(x,y)\;dy\;dx$$
(Note the minus sign on the lower limit of integration.)
You can simplify this a bit by calculating the charge just for the upper half of the disk, and multiplying that by 2. This works because the charge density function is symmetric across the $x$-axis. (That is, $\sigma(x, -y) = \sigma(x, y)$.) Then the integral is $$\color{blue}{2}\int_{x=-1}^{x=1}\int_{\color{blue}{y=0}}^{y=\sqrt{1-x^2}} \sigma(x,y)\;dy\;dx$$
Or similarly you could cut the region into fourths:
$$\color{blue}{4}\int_{\color{blue}{x=0}}^{x=1}\int_{{y=0}}^{y=\sqrt{1-x^2}} \sigma(x,y)\;dy\;dx$$
But probably the best thing do to is to transform the problem into polar coordinates, because it is circularly symmetric, and then the integral becomes $$\int_{\theta=0}^{\theta=2\pi}\int_{r=0}^{r=1} \left(18+r^2\right)\;r\;dr\;d\theta.$$
Starting with your expression (which after cancelling factors can be written)
$$dE = \frac{\sigma z}{2 \epsilon_0} \frac{ada}{(z^2+a^2)^{3/2}}$$
Now integrate to get
$$E = \frac{\sigma z}{2 \epsilon_0} \int_0^R\frac{ada}{(z^2+ a^2)^{3/2}}$$
To solve the integral let $w = z^2 + a^2$ then
$$E = \frac{\sigma z}{4 \epsilon_0} \int_{z^2}^{z^2+R^2}\frac{dw}{w^{3/2}} = \frac{\sigma z}{2 \epsilon_0} \left[-\frac{2}{\sqrt{\omega}}\right]_{z^2}^{z^2+R^2}$$
which gives the correct result
$$E = \frac{\sigma}{2 \epsilon_0} \left(1 - \frac{z}{\sqrt{z^2 + R^2}}\right)$$
Best Answer
$$\frac1{|x|}-\frac1{\sqrt{x^2+R^2}}=\frac{\sqrt{x^2+R^2}-|x|}{|x|\sqrt{x^2+R^2}}=\frac{x^2+R^2-x^2}{|x|\sqrt{x^2+R^2}(\sqrt{x^2+R^2}+|x|)}.$$
For $|x|\gg R$, $\sqrt{x^2+R^2}\approx|x|$ and the expression simplifies to $$\frac{R^2}{2|x|^3}.$$