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

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

def Phillips1(x,t):
    Y=x[0]
    G=x[1]
    return [-a*l*Y+a*G-a*u,-b*G-b*fp*Y]

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

for fp in [0,0.5,2]:
    sol = odeint(Phillips1,[0,0],time)
    x = sol[:,0]
    plt.plot(time,x,label='fp='+str(fp))
plt.legend()
plt.show()