library("deSolve")

A<-c(1,1,1,1)
B<-t(matrix(c(-0.00055, 0.00000,-0.00035,-0.00035,     
               0.00000,-0.00045,-0.00070,-0.00035,     
               0.00050, 0.00200,-0.00400,-0.00150,     
               0.00100, 0.00100, 0.00100,-0.00800),4,4))

podminky<-c(N1=5,N2=5,N3=3,N4=1)
casy<-c(0:200)
soustava<-function(t,prom,param) {
  return(list(c(diag(prom)%*%(param[,1]+param[,2:ncol(param)]%*%prom))))
}
parametry<-cbind(A,B)
metoda<-"lsoda"
barvy<-c("green","blue","red","brown")

reseni<-ode(podminky,casy,soustava,parametry,metoda)

plot(reseni,type="l",col="green",lwd=3,which="N1",ylim=c(0,max(c(reseni[,2:5]),na.rm=TRUE)),xaxt="n",main="Velikost populace",xlab="čas")
for(i in 2:5) {
  lines(reseni[,1],reseni[,i],col=barvy[i-1],lty=1,lwd=3)
}
axis(1,labels=c(0:5)*20,at=c(0:5)*40)