model <- function(t, state, parms) {
  with(as.list(c(state,parms)), {
    dB <- s - d*B
    return(list(dB))  
  }) 
}  

donate <- function(t, y, parms) {
  y["B"] <- y["B"] - parms["Delta"]
  return(y)
}

p <- c(s=1,d=1,Delta=1)
s <- c(B=1) 

run(20,tstep=20/100,events=list(func=donate,time=5),ylab="Red blood cell density",main="(a)",legend=FALSE)
p["Delta"] <- -1
run(20,tstep=20/100,events=list(func=donate,time=5),ylab="Red blood cell density",main="(b)",legend=FALSE)

# Make figure
size <- 5 
pdf("donate.pdf",width=2*size,height=size)
par(par(mar=c(2.6,2.6,1.6,0.2),mgp=c(1.5,0.5,0)),mfrow=c(1,2))
p["Delta"] <- 1
run(20,tstep=20/100,ymax=2,events=list(func=donate,time=5),ylab="Red blood cell density",main="(a)",legend=FALSE)
p["Delta"] <- -1
run(20,tstep=20/100,ymax=2,events=list(func=donate,time=5),ylab="Red blood cell density",main="(b)",legend=FALSE)
dev.off()