model <- function(t, state, parms) {
  with(as.list(c(state,parms)), {
    dB <- m - d*B*(1+(e*B)^n)
    return(list(dB))  
  }) 
}  

p <- c(m=1,d=0.03,e=0.01,n=2)
s <- c(B=35)

run()
continue(newton(run()),x="m",y="B",xmin=0.1,xmax=5,ymax=100)