Following the approach in Dormann et al. (2007), you could do something like this:
N <- 3000
p <- 1/N
# generate some points
set.seed(1234)
x.coord <- runif(N,0,100)
y.coord <- runif(N,0,100)
points <- cbind(x.coord,y.coord)
# distance matrix between points
Dd <- as.matrix(dist(points))
# weights matrix
w <- exp(-p * Dd)
Ww <- chol(w)
# errors
z <- t(Ww) %*% rnorm(N,0,1)
# plot
df <- data.frame(x = x.coord, y = y.coord, z = z)
require(ggplot2)
ggplot(df, aes(x = x, y = y, col = z)) +
geom_point() +
scale_colour_gradient(low="red", high="white")
where variable p controls the size of auto-correlation (here I set it to 1/3000 = 0.000333). p = 0 would be no correlation.
Reference: Dormann, C. F., McPherson, J. M., Araujo, M. B., Bivand, R., Bolliger, J., Carl, G., … Wilson, R. (2007). Methods to account for spatial autocorrelation in the analysis of species distributional data: a review. Ecography, 30(5), 609–628.