I managed to get it to work (I think). Code below:

library(deSolve);
Izhikevich <- function(time, init, parms) {
  with(as.list(c(init, parms)),{
    dv <- (0.04*v^2)+(5*v)+140-u+I;
    du <- a*(b*v-u);
    #if (v>=30) v<-c else v<-u+d;
    list( c(dv, du))
  })}
parms=c( a=0.02, b=0.2, c=-65, d=2, I=4);
times=seq(from=1, to=500, by=0.1);
init=c(v=-65, u=0.2);

root <- function(time, init, parms) init[1]

event <- function(time, init, parms) {
  with(as.list(c(init, parms)), {
  init[2] <- init[1] + d
  init[1] <- c
  return(init)
  })
}

out<-ode(y=init, times=times, func=Izhikevich, parms=parms,
events=list(func=event, root=TRUE), rootfun=root)
plot(out)


Thanks everybody for their help.

Jannetta

	[[alternative HTML version deleted]]

