import matplotlib.pyplot as plt
import numpy as np
from scipy.integrate import odeint

a=4
b=2
l=.25
u=1

def Phillips3(x,t):
    Y=x[0]
    G=x[1]
    Z=x[2]
    return [-a*l*Y+a*G-a*u,Z,-b*Z-b*fi*Y]

time = np.arange(0, 20, 0.01)

legenda=[]
for fi in [0,0.2,0.5,1]:
    sol = odeint(Phillips3,[0,0.00001,0],time)
    x = sol[:,0]
    plt.plot(time,x)
    legenda.append(r'$f_i = {fi}$'.format(fi=fi))
    print((a*l+b)*l-fi)
    

plt.title(r'$\ell(\alpha\ell+\beta) = {g}$'.format(g=l*(a*l+b)))
plt.legend(legenda)
plt.show()