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

plane(xmin=-0.25,xmax=1.5,ymax=2.5,portrait=TRUE,tmax=5,tstep=0.1,grid=7,rtol=1e-12)