model <- function(t, state, parms) {
  with(as.list(c(state,parms)), {
    dR <- R*(1 - R/K) - c1*N*R/(1+b1*R)
    dN <- -aN*N + c1*N*R/(1+b1*R) - c2*M*N/(1+b2*N)
    dM <- -aT*M + c2*M*N/(1+b2*N)
    return(list(c(dR, dN, dM)))  
  }) 
}  

p <- c(b1=6,b2=2,c1=5,c2=0.1,aN=0.4,aT=0.01,K=1)
s <- c(R=1,N=0.1,M=0.01)
plane()
plane(show=names(s),zero=F)
newton(c(R=0.1,N=0.23,M=0),plot=T)
run(traject=T)

# Make a 3D plot:
#install.packages("plot3D")
library(plot3D)
data <- run(1e4,tstep=0.5,table=T)
data <- data[5000:1e4,]    # Take the last 5000 time steps
lines3D(x=data$R,y=data$N,z=data$M,xlab="R",ylab="N",zlab="T",col="red")
# Make a Takens reconstruction:
plot(data$R[1:4999],data$R[2:5000],pch=".")

# Add noise
noiseP="parms[\"K\"]<-abs(rnorm(1,1,0.1))"
data <- run(1e4,tstep=0.5,table=T,after=noiseP)
# Take the last 5000 time steps again and plot

noiseS="if(runif(1)<0.1)state[3]=state[3]+abs(rnorm(1,0,0.1))"
data <- run(1e4,tstep=0.5,table=T,after=noiseS)
# Take the last 5000 time steps again and plot

#install.packages("rgl")
library(rgl)
plot3d(x=data$R,y=data$N,z=data$M,xlab="R",ylab="N",zlab="M",type="l")