t=3/4

# Given equations
f(x,y) = 0.10*x-(.1/10000)*x**2-(.1/10000)*t*x*y # Hardwood
g(x,y) = .25*y-(.25/6000)*y**2-(.25/6000)*t*x*y # Softwood

# Streamline plot
stream=streamline_plot((f(x,y),g(x,y)),(x,0,15000),(y,0,10000))

# Finding where the function quals 0
equil=solve([f(x,y)==0,g(x,y)==0],x,y)
eqpts=[[round(i.rhs(),2) for i in eq] for eq in equil]

# Getting rid of not possible value
eqpts.pop()

# Points on graph
eqplot=list_plot(eqpts,size=30,color='red')

display(stream+eqplot)
table(eqpts)

# Loop through each equilibrium point
for eqpt in eqpts:
    x_val = eqpt[0]
    y_val = eqpt[1]
    
    # Create a zoomed-in streamplot around the equilibrium point
    streams_zoom = streamline_plot((f(x, y), g(x, y)), (x, x_val - 1000, x_val + 1000), (y, y_val - 1000, y_val + 1000), color='gray')
    
    # Put the equilibrium point on the graph
    eqplot_zoom = point(eqpt, size=30, color='purple')
    
    # Display the streamplot and equilibrium point for this eqpt
    display(streams_zoom+eqplot_zoom)