pp_shape = function (r1,K,c,a,Th, C) {
#This is a function plotting two lines of the prey dynamics (growth and predation under given predator level)

#Prey:
#r1 - intrinsic growth rate 
#K - prey carrying capacity

#Predator:
#c - prey conversion rate
#a - encounter rate
#Th - handling time
#C - predator population size

#nostep - number of steps for simulation

nobin=100;

x=NULL;
y1=NULL;
y2=NULL;

for (tt in seq(0,K,by=K/nobin)) {
   x=c(x, tt)
   y1= c(y1, r1*tt*(K-tt)/K)
   y2= c(y2, C*(a*tt^2/(1+a*Th*tt^2)))
}

matplot(x,cbind(y1,y2),xlab='Prey density',ylab='Prey density change',type="l",
        col=c("blue","green"),lty=1,lwd=2,
        main=paste("Input predator density =",as.character(C)))
legend("topleft",lty=1,lwd=2,col=c("blue","green"),legend=c("Prey growth","Predation"))

}