Simulating values from a random variable that is a sum of other random variables $X$ is $\mathcal N(0,4)$, $Y$ is $\mathcal N(0,5)$, $Z = X + Y$
I need to simulate 1000 values for each of these variables, $X$,$Y$,$Z$.
I have simulated 1000 values for both $X$ and 1000 values for $Y$.
When simulating 1000 values for $Z$, should I use the values already simulated for $X$ and $Y$?
Or should I simulate new values for $X$ and $Y$?
 A: I feel like the other answers so far have not been completely clear on why you must reuse the same samples that you obtained for $X$ and $Y$: this is necessary to obtain a sample of $(X,Y,Z)$ which has the correct joint distribution.
If $X$ and $Y$ are independent (that wasn't explicit in the question), then:
 $$\text{Cov}(X,Z) = \text{Cov}(X,X+Y)=\text{Cov}(X,X)=\text{Var}(X)=4$$
We then have $\text{Cor}(X,Z) = 2/3$, so $X$ and $Z$ are definitely not independent. Then, intuitively, the method you use to sample $Z$ must use your existing samples for $X$ and $Y$ in one way or another.
If we compare both approaches (using the same $X$ and $Y$, or new $X$ and $Y$):
set.seed(123)

x <- rnorm(1000, 0, sqrt(4))
y <- rnorm(1000, 0, sqrt(5))
z <- x + y

x_new <- rnorm(1000, 0, sqrt(4))
y_new <- rnorm(1000, 0, sqrt(5))
z_new <- x_new + y_new

par(mfrow=c(1,2))
plot(x,z, main = paste0("Same sample: sample correlation = ", format(cor(x,z),digits=3)))
plot(x,z_new, main = paste0("New sample: sample correlation = ", format(cor(x,z_new),digits=3)))
par(mfrow=c(1,1))

We get the following plots:
On the left, reusing the same sample of $X$ and $Y$ yields a correlation between $X$ and $Z$ which is roughly what we expected (2/3). It's not shown here but it also has the correct joint distribution of $Y$ and $Z$, and the full $(X,Y,Z)$. 
On the right, the correlation is roughly zero. Using new samples of $X$ and $Y$ completely destroys the dependence structure with $Z$.
A: If $Z=X+Y$, and you have $X$ and $Y$, then just sum them, that's the simulation for $Z$.
