This is wrong. A wide region (width $\gg\lambda$) of many, evenly distributed in-phase spherical wave emitters, like the Huygens equivalent of a laser beam's cross section, behaves like a phased antenna array, with the result that destructive interference between the emitters cancels radiation deviating significantly from the beam and reinforces radiation along the beam. So you do indeed get arbitrarily narrow beams from Huygens's theory.
What you might be thinking is that this theory does not cope well with a unidirectionally propagating beam. If we replace the laser beam's cross section in-phase with spherical emitters, a narrow beam is radiated both forwards and backwards (a phased antenna array does this, for example). This is where one must call on an obliquity factor to restore unidirectionality. Fresnel-Huygens theory multiplies the spherical radiation pattern by $\frac{1}{2}(1+\cos\theta)$, where $\theta$ is the angle between the direction of propagation and the line joining the centre of the Huygens wavefront and the point in question.
If you look at the Huygens-Fresnel diffraction integral in the Wikipedia link above, then consider the beam cross-section $B$ in the $x-y$ plane and then also a point $P$ with co-ordinates $(x, y, z)$ in the farfield (i.e. so that $R = \sqrt{z^2+y^2+z^2} \gg w$, where $w$ is the maximum width of the beam cross section, then a Huygens superposition of spherical emitters spread over $B$ will beget a disturbance at $P$ of:
$$\int_B \frac{\exp\left(i\,k\,\sqrt{\left(x-x^\prime\right)^2+\left(y-y^\prime\right)^2+z^2}\right)}{\sqrt{\left(x-x^\prime\right)^2+\left(y-y^\prime\right)^2+z^2}}\,h\left(x^\prime,y^\prime\right)\,\mathrm{d}x^\prime\mathrm{d}y^\prime\approx $$
$$\frac{e^{i\,k\,R}}{R} \int_B \exp\left(-i \frac{k}{R}\left(x\,x^\prime+y\,y^\prime\right)\right)\,h\left(x^\prime,y^\prime\right)\,\mathrm{d}x^\prime\mathrm{d}y^\prime$$
where $h\left(x^\prime,y^\prime\right)$ represents any beam apodisation and aberration (intensity variation and phase, respectively) and I have simply pulled the denominator out of the integral: this is justified because, as a proportion of itself, $\sqrt{\left(x-x^\prime\right)^2+\left(y-y^\prime\right)^2+z^2}$ does not vary much for $\left(x^\prime,y^\prime\right)\in B$ as $R\rightarrow\infty$, but, as a number of wavelengths it varies a great deal, so we must keep the numerator. Thus we see two results:
- The expression is symmetric in $z$, i.e. the same disturbance propagates backwards as it does forwards, so that the beam is not unidirectionally propagating. This is why Fresnel introduced his obliquity factor $K(\chi)$ to remove the backwards propagating wave;
- The farfield beam cross section is given by the Fourier transform of the beam cross section at $B$, but with a scale factor $\frac{k}{R}$ in the Fourier transform's argument, i.e. once we get a reasonable distance from the cross section at $z=0$, the "shape" of the projected beam does not change - it is always equal to a Fourier transform of the beam at $z=0$, but the same shape is dilated by a factor proportional to $R$.
This diffraction scheme is called Fraunhoffer diffraction. So the diffracted beam can be made tightly focussed in the farfield if the cross section at $z=0$ is wide. By dint of the Fourier transform, there is an inverse proportionality between the width of the beam at $z=0$ and that in the farfield. This is exactly the behavior of a high quality laser beam - you will always get a dilation proportional to $R$ - the only way to avoid this is to have a plane wave of infinite extent.
To get some intuition for the result, witness that the above theory gives, for the diffraction of the Gaussian beam defined by:
$$h\left(x^\prime,y^\prime\right) = \exp\left(-\frac{{x^\prime}^2+{y^\prime}^2}{2\sigma^2}\right)$$
the diffracted result at the point $P=(x,y,z)$:
$$\frac{e^{i\,k\,R}}{R} \int_B \exp\left(-i \frac{k}{R}\left(x\,x^\prime+y\,y^\prime\right)\right)\,h\left(x^\prime,y^\prime\right)\,\mathrm{d}x^\prime\mathrm{d}y^\prime = \frac{\sigma^2}{R} \exp\left(-\frac{k^2 \sigma ^2 \left(x^2+y^2\right)}{2 R^2}\right)$$
The link you gave set me searching for a more detailed answer, and I learnt an interesting fact:
Huygens' construction works in 1 and 3 dimensions, BUT NOT IN TWO!
The theory behind this is was first derived by Fresnel and later by Kirchoff - the math is explained in detail in this article. It all comes down to the fact that the wave equation for waves propagating from a point can be written as
$$\frac{\partial^2 \phi}{\partial r^2}-\frac{(n-1)(n-3)}{4r^2}\phi = \frac{\partial^2\phi}{\partial t^2}$$
where $n$ is the dimensionality. For $n=1$ or $n=3$ the second term on the left vanishes, and the expression becomes like that of a one-dimensional wave propagating outwards. For $n=2$, like for the ripples on a pond, there is in fact another term that results in waves "traveling backwards". In that case, the usual Huygens construction does not (quite) work.
Quite subtle stuff - and the fact that the mathematical treatment happened more than 100 years after Huygens' initial publication of his ideas (1679: publication of "Traité de la Lumière"; Fresnel was born in 1788) is something that I had not previously appreciated. The casual treatment it is given in French's book makes sense in the context of a big book covering a lot of ground quickly - but I appreciate you brought this question to my attention; it's more interesting than it looked at first sight.
Best Answer
If you integrate over both $x$ and $y$ you need to take phases into account, otherwise you should get 0 because of destructive interference of circular waves originating half a wavelength away from each other in $y$-direction. If you do that, I assume you should be able to correctly calculate the oscillation of some arbitrary point in the plane, but you do not gain any information, because that oscillation is already given by the plane wave itself, so I do not see why anybody would want to do that.
Another problem is, that when representing the waves by complex exponential functions, if I'm not mistaken, the integral with the phases for the oszillation at the origin would look something like $$ \int dx \int dy ~ \exp \left (i \left(k \left( \sqrt{x^2+y^2} + \underbrace{y}_{\text{phase}}\right) \right) -\omega t\right)~, $$ which seems to be hard to solve.