from sympy import *
import numpy as np


c, t = symbols('c, t')

expr = c*(4 - c)
# Solve c(4 - c) = 0:
equil = solve(expr, c)
print(equil)
ys = np.arange(-2, 7) # how to generate an appropriate array includes equil automatically?
dys = np.array([expr.subs(c, yi) for yi in ys])
print(ys, dys)


# Solve the differential equation y' = y(4 - y):
y = Function('y')
dsols = dsolve(Eq(y(t).diff(t), y(t)*(4 - y(t))), y(t))