library(deSolve) parameters<-c(a=-8/3,b=-10,c=28) state<-c(X=1,Y=1,Z=1) Lorenz <- function(t,state,parameters){ with(as.list(c(state,parameters)),{ # rate of change dX <- a*X + Y*Z dY <-b*(Y-Z) dZ <- -X*Y + c*Y -Z # return the rate of change return(list(c(dX,dY,dZ))) }) } times<- seq(0,100,by=0.01) out<-as.data.frame(ode(y=state,times=times,func=Lorenz,parms=parameters)) head(out) par(mfrow=c(2,2), oma=c(0,0,3,0)) plot(times,out$X, type="l", main="X", xlab="time", ylab="-") plot(times,out$Y, type="l", main="Y", xlab="time", ylab="-") plot(out$X, out$Y, pch=".") plot(out$X, out$Z, pch=".") mtext(outer=TRUE,side=3,"Lorenz model", cex=1.5)