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 <- -2

v <- function(x)
  y0+v2/v1*(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=50)
points(x0,y0,pch=16,col="red")

h <- 10^(-6)
fx <- (f(x0+h,y0)-f(x0,y0))/h
fy <- (f(x0,y0+h)-f(x0,y0))/h