library(pacman)
p_load(deSolve)

a <- 0.005
b <- 0.01 
m <- 0.1
#yini <- c(S=1850,L=100)

Lorenz <- function (t, y, parms) {
#with(as.list(c(S=1800,L=100)), {
list(c(y[1]*a, y[2]*b + (y[1]*m)))
#})
}

out <- ode(y=c(S=1850, L=100), times = seq(0,100), func = Lorenz,
parms = NULL)

tail(out, n=10)