Numerical maximization procedure normal distribution

Firstly, I programmed this gaussian density as an R function which takes the parameters mean and variance (mean=4,sd=1) and calculates the density on a given realization x. Then, passing the density function that I have coded to a (numerical) maximizer, I want to get a maximum close to the parameters which have been used to generate the sample. (mean=4,sd=1).

The question was:

Generate 1000 realizations following N(4,1) and then apply a numerical maximization procedure to determine a maximum likelihood estimator for the unknown parameters.

This is for homework and I would like to know whether my solution is correct or not, and how I can improve. In the case it's wrong I am not asking for a complete correction but just some clues. Hence, my problem is to apply this maximization procedure.

library(stats4)
set.seed(123)
n <- 5000
x <- rnorm(n, mean =4, sd = 1)

likelihood <- function(mu, sigma)
{
a = dnorm(x, mu, sigma)
-sum(log(a))
}
mle(minuslogl =likelihood, start = list(mu = 1, sigma = 1), method = "L-BFGSB",lower = c(-Inf, 0), upper = c(Inf, Inf))


Output:

mu=3.9994305, sd=0.9944726