Is your feature request related to a problem? Please describe.
The lambertw method has to take some pains to appropriately handle numerical overflow in an intermediate calculation.

Describe the solution you'd like
An implementation of [1,2], which reformulates the SDE so that it uses a different form of the Lambert W function and does not require calculating such large intermediate values.

Additionally, according to [3] and also my own brief testing, it is also computationally faster than the lambertw method.

Additional context
[1] "On Calculating the Current-Voltage Characteristic of Multi-Diode Models for Organic Solar Cells": https://arxiv.org/abs/1601.02679
[2] "A Robust Approximation to a Lambert-Type Function": https://arxiv.org/abs/1504.01964
[3] "Solar Cells, Lambert W and the LogWright Functions": https://arxiv.org/abs/2307.08099

(I notice also that @markcampanelli is mentioned in the acknowledgements of [2])

Here is a rough implementation of the method, which I observe to be 2x faster than _lambertw_v_from_i while staying within $\pm$ 1e-10 of its output.

def _g(x, n=2): y = np.where(x <= -np.e, x, 1) y = np.where((-np.e < x) & (x < np.e), -np.e + (1 + np.e) * (x + np.e) / (2 * np.e), y) np.log(x, out=y, where=x >= np.e) for _ in range(n): ey = np.exp(y) ey_plus_one = 1 + ey y_ey_x = y + ey - x y = y - 2 * (y_ey_x) * (ey_plus_one) / ((2 * (ey_plus_one)**2 - (y_ey_x)*ey)) return y def _logwright_v_from_i(current, photocurrent, saturation_current, resistance_series, resistance_shunt, nNsVth): output_is_scalar = all(map(np.isscalar, (current, photocurrent, saturation_current, resistance_series, resistance_shunt, nNsVth))) I, IL, I0, Rs, Rsh, a = \ np.broadcast_arrays(current, photocurrent, saturation_current, resistance_series, resistance_shunt, nNsVth) x = np.log(I0 * Rsh / a) + Rsh * (-I + IL + I0) / a V = -I * Rs - a * np.log(I0 * Rsh / a) + a * _g(x) if output_is_scalar: return V.item() else: return V