1941.sagews

var('y')
MS = MatrixSpace(SR,2)
x_prime(x,y) =  x * ( 2 -x -x^2 -y )
y_prime(x,y) = y * ( 16 -2*x-x^2 -y^2)
f(x,y) = ( 2 -x -x^2 -y )
g(x,y) = ( 16 -2*x-x^2 -y^2)
fx(x,y) = derivative(x_prime(x,y),x)
fy(x,y) = derivative(x_prime(x,y),y)
gx(x,y) = derivative(y_prime(x,y),x)
gy(x,y) = derivative(y_prime(x,y),y)

def comm(x,y):
    A = x * fx(x,y) + f(x,y)
    B = x * fy(x,y)
    C = y * gx(x,y)
    D = y * gy(x,y) + g(x,y)
    return MS.matrix([A,B,C,D])

%latex
\section{The equilibrium points}
$$\begin{array}{rcl}
0 &=& F(x,y)\\
0 &=& G(x,y)
\end{array}$$
Solving this system yields:
    \begin{equation}
    (x_\infty , y_\infty) = \sage{solve([x_prime(x,y)==0, y_prime(x,y)==0],x,y)}
    \end{equation}
We only consider the realistic EQ pts.  Those are the ones that lie in the first quadrant.
        \begin{equation}
    (x_\infty , y_\infty) = \left\lbrace(0,0), (1,0),(0,4) \right\rbrace
    \end{equation}

comm(0,0)

[ 2 0] [ 0 16]

comm(1,0)

[-3 -1] [ 0 13]

comm(0,4)

[ -2 0] [ -32 -128]