model1 <- function(t, state, parms) {
  with(as.list(c(state,parms)), {
    fA <- 1 - A/k
    dJ <- s - d1*J - m*J*fA
    dA <- m*J*fA - d2*A
    return(list(c(dJ, dA)))  
  }) 
}  

model2 <- function(t, state, parms) {
  with(as.list(c(state,parms)), {
    fA <- pmax(0,1 - A/k)
    dJ <- s - d1*J - m*J*fA
    dA <- m*J*fA - d2*A
    return(list(c(dJ, dA)))  
  }) 
} 

model <- model1
p <- c(s=1,k=1,d1=1,m=1,d2=0.5)
s <- c(J=1,A=0.1)
plane(xmax=1.5,ymax=1.5)

# Make Figures
size<-5#inch

pdf("seedlinga.pdf",width=size,height=size)
par(mar=c(0.1,0.1,0.1,0.1),xaxt="n",yaxt="n",ann=F,xaxs="i",yaxs="i")
model <- model1
plane(xmax=1.5,ymax=1.5,main="(a)")
newton(c(J=1,A=0.5),plot=TRUE)
dev.off()

pdf("seedlingb.pdf",width=size,height=size)
par(mar=c(0.1,0.1,0.1,0.1),xaxt="n",yaxt="n",ann=F,xaxs="i",yaxs="i")
model <- model2
plane(xmax=1.5,ymax=1.5,main="(b)")
newton(c(J=1,A=0.5),plot=TRUE)
dev.off()