f <- function(x,y)
  (x^2-3*x*y^2-2*y^3+x-4*y)/(x^2+x*y+y^2+1)
h <- 1.e-6
fx <- function(x,y)
 () f(x+h,y)-f(x-h,y))/(2*h)
fy <- function(x,y)
  (f(x,y+h)-f(x,y-h))/(2*h)

fx(-2,3)
fy(3,-1)
v <- c(4,7)
fv <- function(x,y)
  (f(x+h*v[1],y+h*v[2])- f(x,y))/(h*sqrt(v%*%v))
fv <- c(5,-3)