Phase Portrait of Age-dependent Population Model

#Page 92; Problem-7, determination of Phase Potraits for the age-dependent population
import matplotlib.pyplot as plt
import numpy as np
from scipy.integrate import odeint
import pylab as pl

p=1
a=2
b=3
c=4
# 2D system
def dx_dt(x,t):
  return [x[1]+p*(a-b*p), x[1]*(c+(a-b*p))]

# trajectories
ts=np.linspace(0,10,500)
ic=np.linspace(-3,3,6)

for r in ic:
  for s in ic:
    x0=[r,s]
    xs=odeint(dx_dt,x0,ts)
    plt.plot(xs[:,0],xs[:,1],"r-")

# backward trajectories
ts=np.linspace(0,-10,500)
ic=np.linspace(-3,3,6)
for r in ic:
  for s in ic:
    x0=[r,s]
    xs=odeint(dx_dt,x0,ts)
    plt.plot(xs[:,0],xs[:,1],"r-")

#label the axes and everything esle
plt.xlabel("x",fontsize=15)
plt.ylabel("y",fontsize=15)
plt.tick_params(labelsize=15)
plt.xlim(-3,3)
plt.ylim(-3,3)

#vectorfield

X,Y=np.mgrid[-3:3:20j,-3:3:30j]
u=Y+p*(a-b*p)
v=Y*(c+(a-b*p))
pl.quiver(X,Y,u,v,color='b')
plt.show()