If you want to calculate the likelihood that IQ falls within the range $(a,b)$ you have to integrate the probability density function over that region. You can approximate the integral in R using the following sum:
sum(0.001 * dnorm(seq(-1,1,0.001), mean =0, sd =1))
which is your code, with an additional factor of "0.001."
The article is using a step size of 1 instead of 0.001. Of course multiplying by 1 has no effect, so the article writer did not bother to include it explicitly in their code!
If you make the step size smaller, your integral is more accurate, so total probability will be closer to 1 or 100%. But of course, IQ values are always whole numbers in real life! Increments of 0.001 are a mathematical fiction that arise from pretending that the data are perfectly normal. So it's questionable how much you really gain by insisting on this level of mathematical accuracy.
Alternately you can just use the cumulative distribution function (pnorm or qnorm), where the integral has already been done for you. But personally speaking, I view the density function as a more fundamental concept. In most cases I think in terms of the density function, and treat the cumulative distribution as a useful shortcut to avoid calculating integrals.
