library("deSolve")

a1<-2
a2<-2
a3<-2  

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
  
N10<-100
N20<-200
N30<-300

podminky<-c(N1=N10,N2=N20,N3=N30)
casy<-0:150
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="black",lwd=3,which="N1",ylim=c(0,300))
lines(reseni[,1],reseni[,3],col="green",lty=1,lwd=3)
lines(reseni[,1],reseni[,4],col="brown",lty=1,lwd=3)


A<-c(a1,a2,a3)
B<-matrix(data=c(b11,b21,b31,b12,b22,b32,b13,b23,b33),nrow=3,ncol=3)
N0<-c(N10,N20,N30)


podminky<-c(N=N0)
casy<-0:150
soustava<-function(t,prom,param) {
  return(list(c(diag(prom)%*%(param[,1]+param[,-1]%*%prom))))
}
parametry<-cbind(A,B)
metoda<-"lsoda"

reseni<-ode(podminky,casy,soustava,parametry,metoda)
plot(reseni,type="l",col="black",lwd=3,which="N1",ylim=c(0,300))
lines(reseni[,1],reseni[,3],col="green",lty=1,lwd=3)
lines(reseni[,1],reseni[,4],col="brown",lty=1,lwd=3)