library("deSolve")

a1<-1
a2<-1
a3<-1
b11<--0.01
b12<--0.015
b13<--0.003
b21<--0.003
b22<--0.01
b23<--0.015
b31<--0.015
b32<--0.003
b33<--0.01

podminky<-c(N1=30, N2=50, N3=70)
casy<-0:60
soustava<-function(t,prom,param) {
  return(list(c(prom[1]*(param[1]+param[ 4]*prom[1]+param[ 5]*prom[2]+param[ 6]*prom[3]),
                prom[2]*(param[2]+param[ 7]*prom[1]+param[ 8]*prom[2]+param[ 9]*prom[3]),
                prom[3]*(param[3]+param[10]*prom[1]+param[11]*prom[2]+param[12]*prom[3]))))
}
parametry<-c(a1,a2,a3,b11,b12,b13,b21,b22,b23,b31,b32,b33)
metoda<-"lsoda"

reseni<-ode(podminky,casy,soustava,parametry,metoda)
plot(reseni,type="l",col="red",lwd=3,which="N1",main="3 populace")
lines(reseni[,1],reseni[,3],col="blue",lwd=3)
lines(reseni[,1],reseni[,4],col="yellow",lwd=3)

# vnořené cykly
plot(c(0,60),c(0,70),cex=0,xlab="t",ylab="N")
for(i in 0:5) {
  for(j in 0:5) {
    for(k in 0:5) {
      N10<-20*i
      N20<-20*j
      N30<-20*k
      podminky<-c(N1=N10,N2=N20,N3=N30)
      reseni<-ode(podminky,casy,soustava,parametry,metoda)
      lines(reseni[,1],reseni[,2],col="red",lwd=3)
      lines(reseni[,1],reseni[,3],col="blue",lwd=3)
      lines(reseni[,1],reseni[,4],col="yellow",lwd=3)
    }
  }
}