x <- seq(-10, 10, length= 60)
y <- x
f <- function(x,y) { r <- sqrt(x^2+y^2); 10 * sin(r)/r }
z <- outer(x, y, f) # calculates function on the complete x,y-grid
contour(x, y, z)
filled.contour(x, y, z)
image(x, y, z) # resolution should be increased:
op <- par(bg = "white")
persp(x, y, z, theta = 30, phi = 30, expand = 0.5, col = "lightblue")