f <- function(x,y) x^2+y^2; fxy <- expression(x^2+y^2)
fi <- 20; x1 <- -2; x2 <- 2; y1 <- -2; y2 <- 2; z1 <- -1; z2 <- 8
x <- seq(x1, x2, len = 30);  y <- seq(y1, y2, len = 30)
z <- outer(x, y, f)
FIG <- persp(x, y, z, xlim=c(x1,x2), ylim=c(y1,y2), zlim=c(z1,z2), theta = 120,
   phi = fi, expand = 0.5, col = "yellow",ticktype="detailed",border="red")
lines(trans3d(c(x2,0,0,0,0),c(0,0,y2,0,0),c(0,0,0,0,z2),pmat=FIG),col="blue",lty=2)
dev.new()
contour(x,y,z)
grad <- function(x,y) c(eval(D(fxy,"x")),eval(D(fxy,"y")))
fxy11 <- D(D(fxy,"x"),"x"); fxy12 <- D(D(fxy,"x"),"y")
fxy21 <- D(D(fxy,"y"),"x"); fxy22 <- D(D(fxy,"y"),"y")
H <- function(x,y) eval(fxy11)*eval(fxy22)-eval(fxy12)*eval(fxy21)
grad(0,0)
[1] 0 0
H(0,0)
[1] 4
fxy11
[1] 2

dev.new()
f <- function(x,y) x*y; fxy <- expression(x*y)
fi <- 10; x1 <- -2; x2 <- 2; y1 <- -2; y2 <- 2; z1 <- -5; z2 <- 5
x <- seq(x1, x2, len = 30);  y <- seq(y1, y2, len = 30)
z <- outer(x, y, f)
FIG <- persp(x, y, z, xlim=c(x1,x2), ylim=c(y1,y2), zlim=c(z1,z2), theta = 55,
   phi = 20, expand = 0.5, col = "yellow",ticktype="detailed",border="red")
lines(trans3d(c(x2,0,0,0,0),c(0,0,y2,0,0),c(0,0,0,0,z2),pmat=FIG),col="blue")
dev.new()
contour(x,y,z)
grad <- function(x,y) c(eval(D(fxy,"x")),eval(D(fxy,"y")))
fxy11 <- D(D(fxy,"x"),"x"); fxy12 <- D(D(fxy,"x"),"y")
fxy21 <- D(D(fxy,"y"),"x"); fxy22 <- D(D(fxy,"y"),"y")
H <- function(x,y) eval(fxy11)*eval(fxy22)-eval(fxy12)*eval(fxy21)
grad(0,0)
[1] 0 0
H(0,0)
[1] -1

Vedi qui per approfondimenti "matematici"