model <- function(t, state, parms) {
  with(as.list(c(state,parms)), {
    dB <- r*B*(1 - B/K) - k*N*B/(1 + B/h)
    dN <- s - d*N
    return(list(c(dB, dN)))  
  }) 
}  

p <- c(r=1,K=1,k=20,h=0.1,d=1,s=0.1)
with(as.list(c(s,p)),h*(k*N/r-1))
s <- c(B=0,N=0.1)
plane(ymax=0.25)
newton(c(B=0,N=0.1),plot=TRUE)
newton(c(B=0.1,N=0.1),plot=TRUE)
newton(c(B=0.8,N=0.1),plot=TRUE)
run(tmax=10,tstep=0.1,state=c(B=0.14,N=0.1))
run(tmax=10,tstep=0.1,state=c(B=0.13,N=0.1),add=TRUE)
run(tmax=10,tstep=0.1,state=c(B=0.12,N=0.1),add=TRUE)

#With bone marrow

model <- function(t, state, parms) {
  with(as.list(c(state,parms)), {
    dB <- r*B*(1 - B/K) - k*N*B/(1 + B/h)
    dN <- e*P*(a+(1-a)*B/(i+B)) - d*N
    dP <- s - e*P*(a+(1-a)*B/(i+B))
    return(list(c(dB, dN, dP)))  
  }) 
}  

qss <- function(t, state, parms) {
  with(as.list(c(state,parms)), {
    N <- e*P*(a+(1-a)*B/(i+B))/d
    dB <- r*B*(1 - B/K) - k*N*B/(1 + B/h)
    dP <- s - e*P*(a+(1-a)*B/(i+B))
    return(list(c(dB, dP)))  
  }) 
}  

p <- c(r=1,K=1,k=20,h=0.1,d=1,i=0.2,s=0.1,e=1,a=0.1)
s <- c(B=0,N=0.1,P=1)
run()
plane(odes=qss,state=c(B=1,P=1))
newton(c(B=0,N=0.1,P=1),plot=TRUE,y="P")
newton(c(B=0.1,N=0.1,P=0.2),plot=TRUE,y="P")
newton(c(B=0.8,N=0.1,P=0.2),plot=TRUE,y="P")

run(tmax=10,tstep=0.1,state=c(B=0.1,N=0.1,P=1))
run(tmax=10,tstep=0.1,state=c(B=0.5,N=0.1,P=1),add=TRUE)
run(tmax=10,tstep=0.1,state=c(B=0.9,N=0.1,P=1),add=TRUE)

run(tmax=10,tstep=0.1,state=c(B=0.1,N=0.1,P=0.1))
run(tmax=10,tstep=0.1,state=c(B=0.1,N=0.1,P=0.2),add=TRUE)