supply/demand outcomes for redundant society solution to a system second order bernoulli ode and a first order ode

import numpy as np
from scipy.integrate import odeint
import matplotlib.pyplot as plt

# Define the system of differential equations
def system(y, t, r, K, R0, Rmax, a):
    #dydt = r * y[0] * (1 - y/ K) * (R/ Rmax)
    dydt = r * y[0] * (1 - y[0]/ K) * (y[1]/ Rmax)
    #dRdt = -a * y
    dRdt = -a * y[0]
    return [dydt, dRdt]

# Parameter values
r = 0.1
K = 1000
Rmax = 500
a = 0.05
R = 400

# Initial conditions
y0 = 100
R0 = 400

# Time points
t = np.linspace(0, 100, 1000)  # Time from 0 to 100 units

# Solve the system of equations
solution = odeint(system, [y0, R0], t, args=(r, K, R0, Rmax, a))

# Extract population and resource values
population = solution[:, 0]
resources = solution[:, 1]

# Plot the results
plt.figure(figsize=(10, 6))
plt.plot(t, population, label='Population')
plt.plot(t, resources, label='Resources')
plt.xlabel('Time')
plt.ylabel('Population and Resources')
plt.legend()
plt.title('Population-Resource Dynamics')
plt.grid(True)
plt.show()

population.max()

import pandas as pd
res = pd.DataFrame({'tid':t,'res' :resources})
res[(res['res']<1)&(res['res']>-1)]

res[(res['res']<-199)&(res['res']>-201)]

res[res['tid']<17]