How to do a two dimensional margin distribution?

I am ask to sample and then do a marginal distribution $f(x_1,x_5)$ and plot it. I have the following code which works but dnorm is used for one dimension so I was wondering if i need to change it to dmvnorm.

If so, i change mu=mu.marginal, sigma=sigma.marginal, added a y sample, but dmvnorm says error because of non array input. Anye help is appreciated.

Model of multivariable normal:

mu = c(1,2,6,2,-4)
sigma = c(1,2,1,0.5,2)
rho = diag(rep(1,5))
rho[1,2] = rho[2,1] = 0.4
rho[1,3] = rho[3,1] = -0.3
rho[1,4] = rho[4,1] = -0.7
rho[3,5] = rho[5,3] = 0.2
rho[4,5] = rho[5,4] = 0.5
Sigma = rho * (sigma %o% sigma)

my code:

p = c(1,5)
(mu.marginal = mu[p])
(Sigma.marginal = Sigma[p,p])
# p is one-dimensional: use dnorm() to compute marginal pdf
x = seq(-1,6,by=0.01)

fx = dnorm(x,mean=mu.marginal,sd=sqrt(Sigma.marginal))
ggplot(data=data.frame(x=x,y=fx),mapping=aes(x=x,y=y)) + geom_line(col="blue")

It seems to me you were on the right track with mvtnorm and came close to a solution... I'm not sure how you ran into a non-array input error, but here's what I got with using mvtnorm:

set.seed(123)
dat <- mvtnorm::rmvnorm(1e4, mean = mu.marginal, sigma = Sigma.marginal)
dat <- as.data.frame(dat)

ggplot(dat, aes(x = V1, y = V2)) +
    geom_bin2d()

You can see it's fairly spherical, which is what you'd expect since the off-diagonal elements of Sigma.marginal are 0 (which means that x_1 and x_5 are marginally independently normally distributed...)