f <- function(x,y)
  (27*y^2-26*x*y-46*y+19*x^2-9*x-34)/(16*y^2-48*x*y+16*y+36*x^2-24*x+6)

x0 <- -8
y0 <- 7
#l1 <- 2
#l2 <- -1
r <- 1

l <- function(x)
  y0+l2/l1*(x-x0)

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=20)
points(x0,y0,pch=16,col="red")
abline(h=y0,v=x0,col="black")
lines(x,l(x),col="blue")

h <- 1e-6
fx <- (f(x0+h,y0)-f(x0,y0))/h
fy <- (f(x0,y0+h)-f(x0,y0))/h
fv <- (f(x0+h*l1,y0+h*l2)-f(x0,y0))/
  (h*sqrt(l1^2+l2^2))

fx
fy