Yes, the field is infinite, but it is only log divergent near the plate, so that it is hard to see the divergence numerically. You can see this easily by solving the problem of a uniformly charged infinite plate, which is a 2d problem. Here the charges are uniform along the negative real axis, where the 2d space is imagined to be the complex plane.
This problem can be understood as follows: the 2d electrostatic field of a point charge at the origin, written as a map from C to C can be written in complex form as:
$$ E_x + i E_y = \frac{z} { 2\pi |z|^2 }= \frac{1}{ 2\pi \bar{z}}$$
it points radially outwards. This is a pure antiholomorphic function, except at the origin. It's more familiar to deal with holomorphic functions, so conjugate it!
$$ E_x - iE_y = E(z) = {1\over 2\pi z} $$
Now you want to superpose all the charges on the negative z axis. This is a simple integral:
$$ \int_{-\infty}^0 E(z-a) da = - {1\over 2\pi} \log(z) $$
where I threw away an infinite additive constant (you should think of this as calculating the potential difference between the point z=1 and any other point). This is the function with a given fixed cut discontinuity on the negative real axis.
So at the point $r,\theta$, the electric field is
$$ E_x = - {\rho\over 2\pi} \log(r) $$
$$ E_y = {\rho \theta\over 2\pi} $$
Where I have restored the $\rho$. The part in the y-direction is finite, as your intuition says--- the discontinuity is equal to the charge density (this charge density is the cut discontinuity of the electric field analytic function, which is a way of making it obvious that the electric field goes as the log--- the log function as a constant cut discontinuity). The divergent part is in the x direction, and it is only invisible in the bulk disk because when you get close to the surface, you have cancellations from the left and from the right that wash it out.
So the answer is yes, the E field is log divergent, but only the component in the plane of the disk pointing out. The solution of the disk asymptotes to the plane solution in the near disk limit.
The problem is that you're ignoring the angular dependence of your probe point $\mathbf r$, and that is messing with your integral. If your probe point has cylindrical coordinates $(s,\phi,z)$ and your integration variables are $(s',\phi')$, then the distance between the two is
$$\frac 1 {|\mathbf r-\mathbf r'|}=\frac{1}{\sqrt{s^2-2ss'\cos(\phi-\phi')+s'^2+(z-z')^2}}$$
by the cosine rule (draw it!). If you put this into your integral it will no longer vanish.
(A few pointers on the new integral: the new dependence on $\phi'$ is a bit more complicated. The standard practice is to change variables to $\varphi=\phi'-\phi$. This will leave a simpler denominator, and a factor of $\cos(\varphi+\phi)=\cos(\varphi)\cos(\phi)-\sin(\varphi) \sin(\phi)$ on the numerator. One of the two terms will vanish and the other will yield to a change of variables to $u=\cos(\varphi)$.)
Best Answer
First, I think a little intuition could help. If you imagine a disk with a charge $Q$ get smaller, but all the while keeping the charge Q intact, shouldn't it geometrically approach a point with charge Q? So, in theory, we should expect the field to approach that due to a point charge.
Now, intuition aside, let's go to the mathematics. While the factor $\left[ 1-x/(x^2 + R^2)^{1/2} \right]$ does go to zero as $R\to0$, the charge density $\sigma = Q/(\pi R^2)$ goes to infinity. This is the source of your problem, because you end up with an indeterminate form $\infty \times 0$ while calculating the limit, and not 0.
To evaluate the limit, you could use l'Hôpital's rule, after rewriting your formula as an appropriate fraction (and substituting $\sigma$ for its expression in terms of $R$).
As a bonus, a quick way to do this would be to use this handy approximation, which works for small $y$ values: $$ (1+y)^n \approx 1+ny. $$
You can use this if you factor $x^2$ from the square root: $$ \frac{x}{ \left(x^2 + R^2\right)^{1/2} } = \frac{x}{ \lvert x \rvert \left(1 + (R/x)^2\right)^{1/2} } = \frac{x}{\lvert x \rvert}\left(1+\left(R/x\right)^2 \right)^{-1/2}\approx\frac{x}{\lvert x \rvert} \left( 1 - \frac{R^2}{2x^2}\right). $$
Here I used $y=R/x$ and $n=-1/2$. If $R$ goes to 0, then $y$ goes to 0, and this approximation gets better. Replacing in the expression you have given, we end up with $$ \vec E = \frac{Q}{4 \pi \epsilon x^2} \frac{x \hat \imath}{\lvert x \rvert} = \frac{Q}{4 \pi \epsilon x^2} \frac{\vec x}{\lvert x \rvert}, $$ which is the expression for a field due to a point charge. (Notice that the term $\vec x / \lvert x \rvert$ only gives you the direction of the field, but doesn't change its magnitude.)
Edit: if you try to do the calculations for $x<0$ you'll end up in trouble. The actual formula for the electric field should be $$ \vec E = \frac{\sigma}{2\epsilon}\left[ \frac{x}{\lvert x \rvert} - \frac{x}{\left(x^2 + R^2 \right)^{1/2}} \right] \hat \imath, $$ which you can see if you follow the derivation of the equation.