model <- function(t, state, parms) {
  with(as.list(c(state,parms)), {
    dx <- -4*y
    dy <- 4*x - x^2 - 0.5*y
    return(list(c(dx, dy)))  
  }) 
}  
s <- c(x=0,y=0)
p <- NULL
plane(xmin=-1,xmax=5,ymin=-1,ymax=10,vector=TRUE)
newton(c(x=0,y=0),jacobian=TRUE,vector=TRUE,plot=TRUE)
newton(c(x=4,y=0),jacobian=TRUE,vector=TRUE,plot=TRUE)

model <- function(t, state, parms) {
  with(as.list(c(state,parms)), {
    dx <- 9*x + y^2
    dy <- x - y
    return(list(c(dx, dy)))  
  }) 
}
s <- c(x=0,y=0)
p <- NULL
plane(xmin=-10,xmax=1,ymin=-10,ymax=1,vector=TRUE)
newton(c(x=0,y=0),jacobian=TRUE,vector=TRUE,plot=TRUE,atol=1e-20)
newton(c(x=-9,y=-9),jacobian=TRUE,vector=TRUE,plot=TRUE)

model <- function(t, state, parms) {
  with(as.list(c(state,parms)), {
    dx <- 2*x - x*y
    dy <- -y + y^2*x
    return(list(c(dx, dy)))  
  }) 
}
s <- c(x=0,y=0)
p <- NULL
plane(xmin=-1,xmax=3,ymin=-1,ymax=3,vector=TRUE)
newton(c(x=0,y=0),jacobian=TRUE,vector=TRUE,plot=TRUE,atol=1e-20)
newton(c(x=0.5,y=2),jacobian=TRUE,vector=TRUE,plot=TRUE)