I have this system of ODE:
x˙ = x(3 − 2x − y)
y˙ = y(2 − x − y)
and so far I managed to draw the phase portrait

def phase():
#fixed points
fixed_points = np.array([[0., 0.], [0., 2.], [1.5, 0.], [1., 1.]])
x = np.linspace(0, 3, 10)
y = np.linspace(0, 3, 10)
X, Y = np.meshgrid(x, y)
#function
Xdot = X * (3 - 2 * X - Y)
Ydot = Y * (2 - X - Y)
fig, ax = plt.subplots()
plt.title("Phase portrait")
plt.ylabel("y")
plt.xlabel("x")
#X and Y axix
x1 = np.linspace(0, 3)
y2 = np.linspace(0, 3)
x2 = y2 * 0 # x=0
y4 = x1 * 0 # y=0
plt.plot(x1, y4, color="k")
plt.plot(x2, y2, color="k")
ax.streamplot(X, Y, Xdot, Ydot, density = 0.8)
ax.scatter(*fixed_points.T, color="r")
fig.show()
What I'm asked to do now is to draw the stable and unstable manifolds of the fixed points. I managed to calculate them analytically, but I really don't know how to draw them numerically. From what I know they should follow somehow the lines in the phase portrait, but I'm really struggling trying to sketch them.