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

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

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

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

time = np.arange(0, 10, 0.01)
colors=['black','blue','red']
i=0
for fp in [0,0.5,2]:
    sol = odeint(Phillips1,[0,0],time)
    x = sol[:,0]
    plt.plot(time,x,color=colors[i],label='fp='+str(fp))
    i+=1
    
i=0
for fp in [0,0.5,2]:
    sol = odeint(Phillips2,[0,0],time)
    x = sol[:,0]
    plt.plot(time,x,'--',color=colors[i],label='fp='+str(fp)+' fd='+str(fd))
    i+=1
plt.legend()
plt.show()