is there an R function for finding the local minima of a bivariate function?

I have the following function:

I am interested in finding all the 4 local minima of this bivariate function using code in R. How can I go about it?

If you are interested in numerical optimization, you have several approaches possible. The most direct one is to use optim. By default, this is a Nelder-Mead simplex method but others are implemented.

You will need to start from different starting values to converge to different end points. I can propose you the following:

func <- function(a){
  x <- a[1]
  y <- a[2]
  return(
     0.5*(x^4 - 16*x^2 + 5*x + y^4 - 16*y^2 + 5*y)
  )
}

t0 <- rnorm(100, sd = 20)
t1 <- rnorm(100, sd = 20)

points <- do.call(rbind, lapply(1:100, function(i) optim(par = c(t0[i],t1[i]), fn = func)$par))

And if you want to see graphically your solutions:

library(ggplot2)
ggplot(data.frame(points)) + geom_point(aes(x = X1, y = X2))

You have four local minima in this output