library(deSolve)
parameters <- c(s = 10, r = 25, b = 8/3)
state <- c(X = 0, Y = 1, Z = 1)

Lorenz <- function(t, state, parameters) {
  with(as.list(c(state, parameters)), {
    dX <- s * (Y - X)
    dY <- X * (r - Z) - Y
    dZ <- X * Y - b * Z
    list(c(dX, dY, dZ))
  })
}

times <- seq(0, 50, by = 0.01)
out <- ode(y = state, times = times, func = Lorenz, parms = parameters)
par(oma = c(0, 0, 3, 0))
plot(out, xlab = "time", ylab = "-")
plot(out[, "Y"], out[, "Z"], pch = ".", type = "l")
mtext(outer = TRUE, side = 3, "Lorenz model", cex = 1.5)

str(out)

##  deSolve [1:5001, 1:4] 0 0.01 0.02 0.03 0.04 0.05 0.06 0.07 0.08 0.09 ...
##  - attr(*, "dimnames")=List of 2
##   ..$ : NULL
##   ..$ : chr [1:4] "time" "X" "Y" "Z"
##  - attr(*, "istate")= int [1:21] 2 5462 10914 NA 5 6 0 68 23 NA ...
##  - attr(*, "rstate")= num [1:5] 0.01 0.01 50 0 44.72
##  - attr(*, "lengthvar")= int 3
##  - attr(*, "class")= chr [1:2] "deSolve" "matrix"
##  - attr(*, "type")= chr "lsoda"

dev.off()

## null device 
##           1

library(scatterplot3d)
scatterplot3d(x=out[,2],y=out[,3],z=out[,4],grid=T,
              type="l",color="blue", #angle=80,
              box=FALSE,highlight.3d=F,
              xlab="x", ylab="y", zlab="z",
              main="Rzut 3d atraktora Lorenza\ns = 10, r = 15, b = 8/3, t_0 = 0, t_end=50")