How do you fit an implicit function on Python?

This is the equation I have to fit. The function fit can be done on Origin but the process is tedious when working with several datasets of V_bias and I_out. I'm still a newbie, but I believe it can be much faster when it is done in Python.

What comes to mind is scipy.odr or curve_fit, but this is my first encounter with an implicit function, thus I'm really struggling with this.

def func(var, params):
    x, y = var
    a, b, c, d = params
    return a - b * (np.exp((q *(x + y * c)) / m * kb * T) - 1) - (x + y * c) / d

Model = scipy.odr.Model(func, implicit=True)
Data = scipy.odr.Data([x,y], 1)
Odr = scipy.odr.ODR(Data, Model, [a,b,c,d], maxit=10000)

ValueError: too many values to unpack (expected 2)

Here's an attempt I borrowed from Implicit fitting with scipy.odr

But it gives me that ValueError. How do I proceed?