I do not really know how much this can help you, but maybe you can use this code for rough sketches of the dynamics of planar vector fields. I know that it does not have the functionality that you may really wish for, but it may come in handy to guide the analysis in the right direction.
import numpy as np
import matplotlib.pyplot as plt
This is not surprising. For one, Broyden is not very stable for highly ill-conditioned nonlinear systems, and having a stiff ODE means that the nonlinear system will be ill-conditioned. This is why stiff ODE solvers don't use Broyden but use Newton iterations (well, a special form of Quasi-Newton which is Newton iteration but sometimes reusing the Jacobian a ...
Since performance is not a primary concern, just use the backward Euler finite difference approximation in time, full Newton method with a direct solver at each time step, and a central difference approximation to the Jacobian.
Since the Jacobian is only a 2x2 matrix, a quasi-Newton solution approach is not useful.