def ode_schnakenberg(t, y, a_prod, b_prod):
    """Derivatives to be called into solve_ivp
    
    This returns an array of derivatives y' = [A', B'], for a given
    state [A, B] at a time t. This is based on the classical 
    Schnakenberg system.

    Params:
    t [float] - the time at which the derivative is evaluated
    y [array] - the current state of the system, in the form [A, B]
    """

    return [y[0]**2 * y[1] - y[0] + a_prod, -y[0]**2 * y[1] + b_prod]