lv2_dist_plot = function(r1,k1,n01,r2,k2,n02,alfa12,alfa21,frequ,int1,int2,nostep,animation) {
#This function plots dynamics of LV system with disturbance and isolines together

#r1, r2 - intrinsic growth rates for species 1 and 2
#k1, k2 - carrying capacities
#n01, n02 - initial population sizes
#alfa12, alfa21 - competition coefficients
#frequ - how often to disturb (no of time steps between disturbance events)
#int1, int2 - proportion of population of species 1 and 2 left after disturbance  
#nostep - number of steps for simulation
#animation - whether to animate

  x=lv2_dist(r1,k1,n01,r2,k2,n02,alfa12,alfa21,frequ,int1,int2,nostep);
  y=lv2_lines(k1,k2,alfa12,alfa21);

  matplot(y[,1],y[,2:3], xlim=c(0, max(c(y[,1],x[,1]),na.rm=T)), ylim=c(0, max(c(y[,2:3],x[,2]),na.rm=T)),
          type="l", col=c("blue","green"),lwd=2,lty=1,
          xlab="Population size of the species 1",ylab = "Population size of the species 2")
  if (animation > 0) {
    for (i in 1:(nostep-1)) { 
      lines(c(x[i,1],x[i+1,1]),c(x[i,2],x[i+1,2]),col="red",lwd=2);
      Sys.sleep(0.1);
    }
  } 
  else {
    lines(x[,1],x[,2],col="red",lwd=2);
  }
}