def F(X):
    x,y = X
    return np.array( (y**2 - 2*x - 3, y-x ))

plt.subplot(111,aspect=1)
fieldplot(F,-3,5,-3,5,color='k',alpha=0.1)
x = np.linspace(-3,5,51)
y = np.linspace(-3,5,51)
X,Y = np.meshgrid(x,y)

plt.grid();

x,y = sp.symbols('x,y')
f,g = F((x,y))
J = sp.Matrix([ [sp.diff(f,x),sp.diff(f,y)],
                [sp.diff(g,x),sp.diff(g,y)]])
display(J)
J1 = J.subs({x:-1,y:-1})
display(J1)
l = sp.symbols('lambda')
display(((J1-l*sp.eye(2))).det())
eig = J1.eigenvects()
display(eig)
(l,m1,v1),(l2,m2,v2) = eig
display(v1[0].evalf())
v = np.array(v1[0].evalf(),dtype=float).reshape(2)
v /= np.linalg.norm(v)
v /=2
display(v)
plt.plot([-1,-1+v[0]],[-1,-1+v[1]],'k',lw=3)
plt.plot([-1,-1-v[0]],[-1,-1-v[1]],'k',alpha=0.25,lw=3)
v = np.array(v2[0].evalf(),dtype=float).reshape(2)
v /= np.linalg.norm(v)
v /=2
display(v)
plt.plot([-1,-1+v[0]],[-1,-1+v[1]],'k',lw=3)
plt.plot([-1,-1-v[0]],[-1,-1-v[1]],'k',alpha=0.25,lw=3)

J2 = J.subs({x:3,y:3})
display(J2)
l = sp.symbols('lambda')
display(((J2-l*sp.eye(2))).det())
eig = J2.eigenvects()
display(eig)

cmap = matplotlib.colormaps['rainbow']
y0s = np.linspace(-3,5,17) 
colors = cmap(np.linspace(0,1,len(y0s)))
plt.plot([-1,3],[-1,3],'ro',ms=15,alpha=0.5)
for i,y0 in enumerate(y0s):
    phaseportrait(F,[(0,y0,-1.2,1.2)], alpha=0.95,color=colors[i],lw=2.5)
for i,y0 in enumerate(y0s):
    phaseportrait(F,[(-2,y0,-1.2,1.2)],alpha=0.95,color=colors[i],lw=2.5)#,color='m')
for i,y0 in enumerate(y0s):
    phaseportrait(F,[(3,y0,-1.2,1.2)], alpha=0.95,color=colors[i],lw=2.5)#,color='m')
plt.xlim(-3,5)
plt.ylim(-3,5)

plt.contour(X,Y,F((X,Y))[0],levels=[0])
plt.contour(X,Y,F((X,Y))[1],levels=[0]);
plt.title('phase portrait')
plt.xlabel('x',fontsize=15)
plt.ylabel('y',fontsize=15,rotation=0);