diff --git a/pvlib/pvsystem.py b/pvlib/pvsystem.py
index a8fc0c7258..9144106e79 100644
--- a/pvlib/pvsystem.py
+++ b/pvlib/pvsystem.py
@@ -2662,7 +2662,7 @@ def v_from_i(current, photocurrent, saturation_current, resistance_series,
         V = _singlediode.bishop88_v_from_i(*args, method=method.lower())
         if all(map(np.isscalar, args)):
             return V
-        shape = np.broadcast_shapes(*map(np.shape, args))
+        shape = _singlediode._shape_of_max_size(*args)
         return np.broadcast_to(V, shape)
 
 
@@ -2744,7 +2744,7 @@ def i_from_v(voltage, photocurrent, saturation_current, resistance_series,
         current = _singlediode.bishop88_i_from_v(*args, method=method.lower())
         if all(map(np.isscalar, args)):
             return current
-        shape = np.broadcast_shapes(*map(np.shape, args))
+        shape = _singlediode._shape_of_max_size(*args)
         return np.broadcast_to(current, shape)
 
 
diff --git a/pvlib/singlediode.py b/pvlib/singlediode.py
index c82a69bb53..a00554f49e 100644
--- a/pvlib/singlediode.py
+++ b/pvlib/singlediode.py
@@ -596,6 +596,11 @@ def fmpp(x, *a):
         return bishop88(vd, *args)
 
 
+def _shape_of_max_size(*args):
+    return max(((np.size(a), np.shape(a)) for a in args),
+               key=lambda t: t[0])[1]
+
+
 def _prepare_newton_inputs(x0, args, method_kwargs):
     """
     Make inputs compatible with Scipy's newton by:
@@ -620,7 +625,7 @@ def _prepare_newton_inputs(x0, args, method_kwargs):
     """
     if not (np.isscalar(x0) and all(map(np.isscalar, args))):
         args = tuple(map(np.asarray, args))
-        x0 = np.broadcast_to(x0, np.broadcast_shapes(*map(np.shape, args)))
+        x0 = np.broadcast_to(x0, _shape_of_max_size(x0, *args))
 
     # set abs tolerance and maxiter from method_kwargs if not provided
     # apply defaults, but giving priority to user-specified values