# -*- coding: utf-8 -*-
from scipy.integrate import odeint
from pylab import * # for plotting commands
import numpy as np

    
def od_deriv(y,t,a,b,d):
    u = y[0]
    v = y[1]
    du = u*(1-u) - a * u*v/(u+d)
    dv = b *v*(1-v/u)
    return array([du, dv])
    
def pole_wektorowe(zakres_x,zakres_y,a,b,d ):
    u,v = meshgrid( zakres_x, zakres_y )
    U = u*(1-u) - a * u*v/(u+d)
    V = b *v*(1-v/u)
    quiver(u,v,U,V)
    
# warunek stabilności
    
figure(3)
d,a = meshgrid(arange(0.1,2,0.02),arange(0.1,3,0.1))
u_s = (1 - d-a +np.sqrt((1-d-a)**2 +4*d))/2
lewa = u_s*(1-2*u_s-d)/(u_s+d)
pcolor(d,a,lewa)
xlabel('d')
ylabel('a')
colorbar()

show()

# całkowanie
time = linspace(0.0,1000.0,100000) 
yinit = array([1.5,.5]) # initial values 
a = 2.
b = 0.01
d= 0.4
for i in range(1):
    y = odeint(od_deriv,yinit+i*0.1,time, args =(a,b,d),atol=1e-13,rtol = 1e-13)
    figure(1) 
    plot(time,y[:,0]) # y[:,0] is the first column of y 
    plot(time,y[:,1])
    xlabel('t') 
    ylabel('y') 
    
    figure(2)
    plot(y[:,0],y[:,1]) # y[:,0] is the first column of y 
    pole_wektorowe(arange(0.1,1,0.05),arange(0.15,0.3,0.01),a,b,d)
show()