# Source grind.R and cube.R first!

model <- function(t, state, parms) {
  with(as.list(c(state,parms)), {
     dN1 <- r*N1*(1-N1-a*N2)
     dN2 <- r*N2*(1-N2-a*N1-a*N3)
     dN3 <- r*N3*(1-N3-a*N2)
    return(list(c(dN1, dN2, dN3)))  
  }) 
}  

p <- c(r=1,a=0.5)
s <- c(N1=1,N2=0.1,N3=1)
cube(xmax=2,ymax=2,zmax=2)
plotdev(theta=50,phi=25)

size<-5#inch

pdf("invasion.pdf",width=size,height=size)
par(mar=c(0,0,0,0))
cube(xmax=2,ymax=2,zmax=2,theta=50,phi=25)
f <- newton(state=c(N1=0,N2=0,N3=0))
points3D(x=f[1],y=f[2],z=f[3],add=TRUE,pch=1)
f <- newton(state=c(N1=1,N2=0,N3=0))
points3D(x=f[1],y=f[2],z=f[3],add=TRUE,pch=1)
f <- newton(state=c(N1=0,N2=1,N3=0))
points3D(x=f[1],y=f[2],z=f[3],add=TRUE,pch=1)
f <- newton(state=c(N1=0,N2=0,N3=1))
points3D(x=f[1],y=f[2],z=f[3],add=TRUE,pch=1)
dev.off()