x0 <- 0
y0 <- 0
r <- 10

f <- function(x,y)
  x^3+5*y^3-3*x*y^2-4*x^2-2*x*y+5*x-2*y

x <- seq(x0-r,x0+r,r/100)
y <-seq(y0-r,y0+r,r/100)
z <- outer(x,y,f)

contour(x,y,z,asp=1,nlevels=50)