from scipy.integrate import odeint
from pylab import * 
    
def N_deriv(y,t,r,K,E):
    N = y[0]
    dN = r*N*(1-N/K) -E*N
    return array([dN])
    
time = linspace(0.0,50.0,1000)

Ninit = array([0.5]) 
r = 2.0
K = 1.
E0 = 0.5 

N_skala = np.arange(0,0.7, 0.01)

for i in range(1,10):
    E = E0 + i*0.15 
    N_ini = odeint(N_deriv,Ninit,time, args =(r,K,E),atol=1e-13,rtol = 1e-13)
    N_dalej = odeint(N_deriv,N_ini[-1] -0.02,time, args =(r,K,E),atol=1e-13,rtol = 1e-13)
    N = np.concatenate((N_ini,N_dalej),axis=0)
    t = np.concatenate((time,time[-1]+time),axis=0)
    #N_bardziej = odeint(N_deriv,N_ini[-1],time, args =(r,K,E+0.1),atol=1e-13,rtol = 1e-13)
    #N_odpuszczam = odeint(N_deriv,N_bardziej[-1],time, args =(r,K,E),atol=1e-13,rtol = 1e-13)
    #N = np.concatenate((N_ini,N_bardziej,N_odpuszczam),axis=0)
    #t = np.concatenate((time,time[-1]+time,2*time[-1]+time),axis=0)
    
    figure(1)
    subplot(1,3,1)
    plot(t,N) 
    xlabel('t') 
    ylabel('N')
    text(t[-1]-20,N[-1], 'E = '+str(E) )
    title('liczebnosc populacji')
    

    subplot(1,3,2)
    zysk  = N*E
    plot(t, zysk)
    text(t[-1]-20,zysk[-1], 'E = '+str(E) )
    ylim([0,0.55])
    xlabel('t')
    title('zysk')
    
    subplot(1,3,3)
    f = r*N_skala*(1-N_skala/K) -E*N_skala
    plot(N_skala, f) 
    ylim([0,0.25])
    title('f(N)')
    xlabel('N')
show()