Visualising OLS estimator sampling distribution

I am trying to understand the idea of the sampling distribution of an OLS estimator. To do so I am trying to visualise this distribution obtained with simulations (I use R for this).

So I generate the universe of observations with 100,000 randomly distributed observations of my independent variable univ.x. Then I also generate the universe of observations for the dependent variable (univ.y) according to the following (exact) linear relationship: univ.y=1.5+3*univ.x.

> univ.x <- rnorm(100000,0,1)
> univ.y <- 1.5 + 3*univ.x
> summary(lm(univ.y~univ.x))

Call:
lm(formula = univ.y ~ univ.x)

Residuals:
Min         1Q     Median         3Q        Max
-2.900e-15 -3.000e-16 -2.000e-16  0.000e+00  8.498e-12

Coefficients:
Estimate Std. Error   t value Pr(>|t|)
(Intercept) 1.50e+00   9.75e-17 1.538e+16   <2e-16 ***
univ.x      3.00e+00   9.73e-17 3.083e+16   <2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 3.083e-14 on 99998 degrees of freedom
Multiple R-squared:      1, Adjusted R-squared:      1
F-statistic: 9.507e+32 on 1 and 99998 DF,  p-value: < 2.2e-16

So far so good. At this point, I want to sample 1,000 observations from the universe and I would expect to find the same linear relationship with some random error. But I don't, in fact I find only the intercept but the slope is 0 and not 3. I have also repeated this 100 times with the same result.

What am I missing here?

Here is the code:

> x <- sample(univ.x, 1000)
> y <- sample(univ.y, 1000)
> summary(lm(y ~ x))

Call:
lm(formula = y ~ x)

Residuals:
Min      1Q  Median      3Q     Max
-9.9094 -2.0161 -0.0021  1.9563  8.9696

Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 1.505428   0.096626  15.580   <2e-16 ***
x           0.009807   0.094319   0.104    0.917
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 3.055 on 998 degrees of freedom
Multiple R-squared:  1.083e-05, Adjusted R-squared:  -0.0009912
F-statistic: 0.01081 on 1 and 998 DF,  p-value: 0.9172

Visualising the sampling distribution of OLS estimates [solved]

Thanks to the comments I solved my problem. I leave the code here for interested users and learners.

univ.x <- rnorm(100000,0,1)
univ.y <- 1.5 + 3*univ.x + rnorm(100000,0,3)
summary(lm(univ.y~univ.x))
plot(univ.x,univ.y)

b <- vector(mode = "numeric", 1000)
for (i in 1:1000){
samp <- sample(1:100000, 100)
x <- univ.x[samp]
y <- univ.y[samp]
model <- lm(y~x)
b[i] <- coef(model)
}
hist(b)

You are sampling 1000 of the $x_i$ at random and 1000 of the $y_i$ at random, so you are not necessarily sampling 1000 pairs $(x_i,y_i)$ that fulfill the linear relationship, or that are connected in any way.
Moreover, the data you generated has no randomness in the $y_i$ respect to the $x_i$. Once you knew one observation of $x_i$ you would know the corresponding $y_i$ for sure, right? This is not the real-life setting that OLS or regression is designed for.
The true regression equation you would get from real life data is $$y_i=\beta x_i+\varepsilon_i,$$ where $\varepsilon_i\sim N(0,\sigma2)$. So to obtain randomness you also need to generate the random errors $\varepsilon_i$, all of them drawn from the same distribution $N(0,\sigma2)$.