import matplotlib.pyplot as plt
import numpy as np

def dydt(y, t):
    return (5-y) * (y+1)**2 * np.sin(y)

# Create mesh points
npts = 20
t, y = np.linspace(0, 1, npts), np.linspace(-7, 7, npts)
T, Y = np.meshgrid(t, y)

# At each mesh point, create a unit vector with correct slope
dy = dydt(Y, T); dt = np.ones(dy.shape)
norm = np.sqrt(T**2 + Y**2) 
dyu, dtu = dy/norm, dt/norm

# Plot
plt.subplots(figsize=(7,5))
plt.streamplot(T, Y, dtu, dyu, linewidth=.5, arrowsize=.7, density=[2])
plt.xlabel('t')
plt.ylabel('y(t)')
plt.yticks(np.arange(-7, 8, 1))
plt.title(r"$Stream\ Plot\ for\ \ y'=(5-y)(y+1)^2\sin(y)$")