Why is the Tukey's IQR not used in the R program? I was writing an R Markdown for calculating measures of spread and was surprised to find that the IQR function does not use the IQR calculation I am used to using, which is:
$$
\text{IQR} = \text{Q3} - \text{Q1}
$$
So with a vector like this:
$$
x = \begin{bmatrix} 10,33,50,341,987,2006,2008 \end{bmatrix}
$$
Q3 should be $2006$ and Q1 should be $33$. Therefore:
$$
\text{IQR} = 2006 - 33 = 1973
$$
However, when running this in R, this does not appear to be how it is derived in the software.
x7 <- c(10,33,50,341,987,2006,2008)
IQR(x7)

Which returns this:
[1] 1455

Looking through the help page, it gives this brief explanation:

Note that this function computes the quartiles using the quantile
function rather than following Tukey's recommendations, i.e., IQR(x) =
quantile(x, 3/4) - quantile(x, 1/4).
For normally N(m,1)N(m,1) distributed XX, the expected value of IQR(X)
is 2*qnorm(3/4) = 1.3490, i.e., for a normal-consistent estimate of
the standard deviation, use IQR(x) / 1.349.

My question is: why would this method be the preferred default if a lot of data is not normally distributed?

Ben Bolker's answer explains why the quantile() function, and thus the IQR, may not return the values you expect. However, I think you've misunderstood the documentation about normality.
Only the "Type 9" method for quantile and IQR is intended for normally-distributed data, where it provides a median-unbiased estimate (docs here). The default is Type 7.
Instead, the IQR documentation is hinting at a method for robustly estimating the standard deviation of samples from a normal distribution. You might be tempted to do so by plugging the values directly into the (Bessel-corrected) definition:
$$
s = \sqrt{\frac{1}{n-1} \sum_{i=0}^{N} \left(x_i - \bar{x}\right)^2}
$$
However, this can go badly wrong with noisy data. The mean $\bar{x}$ is extremely sensitive to outliers. In fact, even in a single outlier is enough to make $\bar{x}$ take on any value. Robust methods can provide more sensible estimates of the SD.
One popular way to do so for normally distributed data is to calculate the IQR and divide it by 1.349. This trick exploits the fact that, for a standard normal distribution, the 1st and 3rd quartiles occur at $\approx \pm 0.67449$). Thus, if you know the data is normal(ish) and the IQR that it spans, you can divide that IQR by 0.67449 * 2, or 1.349, and approximate the SD. There are other robust methods based around the median absolute deviation or pairwise differences. These are useful because the order statistics are much more resilient: changing one value has almost no effect on the median. For example, you'd need to change >50% of them instead to set it to an arbitrary value.
