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
         list(c(dX, dY, dZ))
       }) # end with(as.list ...
};
Span=0.006
SimTime=100;
IterTimes=seq(1,SimTime/Span,by=1);
out=matrix(data = NA, nrow = length(IterTimes), ncol = 4)
for (i in IterTimes){
  times <- seq(0, Span, by = Span);
  out1 <- ode(y = state, times = times, func = Lorenz, parms = parameters);
  nrosOut1=dim(out1)[1];
  state <- c(out1[nrosOut1,2],out1[nrosOut1,3], out1[nrosOut1,4]);
  out[i,]=c(Span*(i-1),state);
}

head(out);
tail(out);
par(oma = c(0, 0, 3, 0),lty="dotdash",col=rgb(0, 0, 1))#Set or Query Graphical Parameters
plot(out[,1], out[,2],xlab = "time", ylab = "x(t)",type="l",col=rgb(0, 0, 1))
mtext(outer = TRUE, side = 3, "Lorenz model", cex = 1.5)
plot(out[,1], out[,3],xlab = "time", ylab = "y(t)",type="l",col=rgb(0, 1, 0))
mtext(outer = TRUE, side = 3, "Lorenz model", cex = 1.5)
plot(out[,1], out[,4],xlab = "time", ylab = "y(t)",type="l",col=rgb(1, 0, 0))
#plot(out[, "3"], out[, "4"], pch = ".")
mtext(outer = TRUE, side = 3, "Lorenz model", cex = 1.5)