f <- function(x, y) { (x^2+y-11)^2 + (x+y^2-7)^2 }

x <- seq(-6, 6, length = 100)
y <- x

z <- outer(x, y, f)

library(plotly)
plot_ly(x = x, y = y, z = ~z) %>% add_surface()