diff --git a/docs/sphinx/source/whatsnew/v0.9.1.rst b/docs/sphinx/source/whatsnew/v0.9.1.rst
index 4ddcb9bd51..80a47ff0b4 100644
--- a/docs/sphinx/source/whatsnew/v0.9.1.rst
+++ b/docs/sphinx/source/whatsnew/v0.9.1.rst
@@ -24,6 +24,10 @@ Bug fixes
   and ``surface_azimuth`` inputs caused an error (:issue:`1127`, :issue:`1332`, :pull:`1361`) 
 * Changed the metadata entry for the wind speed unit to "Wind Speed Units" in
   the PSM3 iotools function (:pull:`1375`)
+* Improved convergence when determining the maximum power point using
+  for :py:func:`pvlib.pvsystem.singlediode` with ``method='lambertw'``. Tolerance
+  is determined for the voltage at the maximum power point, and is improved
+  from 0.01 V to 1e-8 V. (:issue:`1087`, :pull:`1089`)
 
 Testing
 ~~~~~~~
diff --git a/pvlib/ivtools/sdm.py b/pvlib/ivtools/sdm.py
index 04263c34c5..013e63f9a6 100644
--- a/pvlib/ivtools/sdm.py
+++ b/pvlib/ivtools/sdm.py
@@ -979,7 +979,8 @@ def _filter_params(ee, isc, io, rs, rsh):
     negrs = rs < 0.
     badrs = np.logical_or(rs > rsh, np.isnan(rs))
     imagrs = ~(np.isreal(rs))
-    badio = np.logical_or(~(np.isreal(rs)), io <= 0)
+    badio = np.logical_or(np.logical_or(~(np.isreal(rs)), io <= 0),
+                          np.isnan(io))
     goodr = np.logical_and(~badrsh, ~imagrs)
     goodr = np.logical_and(goodr, ~negrs)
     goodr = np.logical_and(goodr, ~badrs)
diff --git a/pvlib/tests/test_pvsystem.py b/pvlib/tests/test_pvsystem.py
index 505016a7ce..1141e490e9 100644
--- a/pvlib/tests/test_pvsystem.py
+++ b/pvlib/tests/test_pvsystem.py
@@ -1313,32 +1313,35 @@ def test_singlediode_array():
                               resistance_series, resistance_shunt, nNsVth,
                               method='lambertw')
 
-    expected = np.array([
-        0.        ,  0.54538398,  1.43273966,  2.36328163,  3.29255606,
-        4.23101358,  5.16177031,  6.09368251,  7.02197553,  7.96846051,
-        8.88220557])
-
-    assert_allclose(sd['i_mp'], expected, atol=0.01)
+    expected_i = np.array([
+        0., 0.54614798740338, 1.435026463529, 2.3621366610078, 3.2953968319952,
+        4.2303869378787, 5.1655276691892, 6.1000269648604, 7.0333996177802,
+        7.9653036915959, 8.8954716265647])
+    expected_v = np.array([
+        0., 7.0966259059555, 7.9961986643428, 8.2222496810656, 8.3255927555753,
+        8.3766915453915, 8.3988872440242, 8.4027948807891, 8.3941399580559,
+        8.3763655188855, 8.3517057522791])
+
+    assert_allclose(sd['i_mp'], expected_i, atol=1e-8)
+    assert_allclose(sd['v_mp'], expected_v, atol=1e-8)
 
     sd = pvsystem.singlediode(photocurrent, saturation_current,
                               resistance_series, resistance_shunt, nNsVth)
-
     expected = pvsystem.i_from_v(resistance_shunt, resistance_series, nNsVth,
                                  sd['v_mp'], saturation_current, photocurrent,
                                  method='lambertw')
-
-    assert_allclose(sd['i_mp'], expected, atol=0.01)
+    assert_allclose(sd['i_mp'], expected, atol=1e-8)
 
 
 def test_singlediode_floats():
-    out = pvsystem.singlediode(7, 6e-7, .1, 20, .5, method='lambertw')
-    expected = {'i_xx': 4.2498,
-                'i_mp': 6.1275,
-                'v_oc': 8.1063,
-                'p_mp': 38.1937,
-                'i_x': 6.7558,
-                'i_sc': 6.9651,
-                'v_mp': 6.2331,
+    out = pvsystem.singlediode(7., 6.e-7, .1, 20., .5, method='lambertw')
+    expected = {'i_xx': 4.264060478,
+                'i_mp': 6.136267360,
+                'v_oc': 8.106300147,
+                'p_mp': 38.19421055,
+                'i_x': 6.7558815684,
+                'i_sc': 6.965172322,
+                'v_mp': 6.224339375,
                 'i': None,
                 'v': None}
     assert isinstance(out, dict)
@@ -1346,23 +1349,26 @@ def test_singlediode_floats():
         if k in ['i', 'v']:
             assert v is None
         else:
-            assert_allclose(v, expected[k], atol=1e-3)
+            assert_allclose(v, expected[k], atol=1e-6)
 
 
 def test_singlediode_floats_ivcurve():
-    out = pvsystem.singlediode(7, 6e-7, .1, 20, .5, ivcurve_pnts=3, method='lambertw')
-    expected = {'i_xx': 4.2498,
-                'i_mp': 6.1275,
-                'v_oc': 8.1063,
-                'p_mp': 38.1937,
-                'i_x': 6.7558,
-                'i_sc': 6.9651,
-                'v_mp': 6.2331,
-                'i': np.array([6.965172e+00, 6.755882e+00, 2.575717e-14]),
-                'v': np.array([0., 4.05315, 8.1063])}
+    out = pvsystem.singlediode(7., 6e-7, .1, 20., .5, ivcurve_pnts=3,
+                               method='lambertw')
+    expected = {'i_xx': 4.264060478,
+                'i_mp': 6.136267360,
+                'v_oc': 8.106300147,
+                'p_mp': 38.19421055,
+                'i_x': 6.7558815684,
+                'i_sc': 6.965172322,
+                'v_mp': 6.224339375,
+                'i': np.array([
+                    6.965172322, 6.755881568, 2.664535259e-14]),
+                'v': np.array([
+                    0., 4.053150073, 8.106300147])}
     assert isinstance(out, dict)
     for k, v in out.items():
-        assert_allclose(v, expected[k], atol=1e-3)
+        assert_allclose(v, expected[k], atol=1e-6)
 
 
 def test_singlediode_series_ivcurve(cec_module_params):
@@ -1383,21 +1389,20 @@ def test_singlediode_series_ivcurve(cec_module_params):
     out = pvsystem.singlediode(IL, I0, Rs, Rsh, nNsVth, ivcurve_pnts=3,
                                method='lambertw')
 
-    expected = OrderedDict([('i_sc', array([0., 3.01054475, 6.00675648])),
-                            ('v_oc', array([0., 9.96886962, 10.29530483])),
-                            ('i_mp', array([0., 2.65191983, 5.28594672])),
-                            ('v_mp', array([0., 8.33392491, 8.4159707])),
-                            ('p_mp', array([0., 22.10090078, 44.48637274])),
-                            ('i_x', array([0., 2.88414114, 5.74622046])),
-                            ('i_xx', array([0., 2.04340914, 3.90007956])),
+    expected = OrderedDict([('i_sc', array([0., 3.01079860, 6.00726296])),
+                            ('v_oc', array([0., 9.96959733, 10.29603253])),
+                            ('i_mp', array([0., 2.656285960, 5.290525645])),
+                            ('v_mp', array([0., 8.321092255, 8.409413795])),
+                            ('p_mp', array([0., 22.10320053, 44.49021934])),
+                            ('i_x', array([0., 2.884132006, 5.746202281])),
+                            ('i_xx', array([0., 2.052691562, 3.909673879])),
                             ('v', array([[0., 0., 0.],
-                                         [0., 4.98443481, 9.96886962],
-                                         [0., 5.14765242, 10.29530483]])),
+                                         [0., 4.984798663, 9.969597327],
+                                         [0., 5.148016266, 10.29603253]])),
                             ('i', array([[0., 0., 0.],
-                                         [3.01079860e+00, 2.88414114e+00,
-                                          3.10862447e-14],
-                                         [6.00726296e+00, 5.74622046e+00,
-                                          0.00000000e+00]]))])
+                                         [3.0107985972, 2.8841320056, 0.],
+                                         [6.0072629615, 5.7462022810, 0.]]))])
+
 
     for k, v in out.items():
         assert_allclose(v, expected[k], atol=1e-2)
@@ -1414,7 +1419,7 @@ def test_singlediode_series_ivcurve(cec_module_params):
                                          method='lambertw').T
 
     for k, v in out.items():
-        assert_allclose(v, expected[k], atol=1e-2)
+        assert_allclose(v, expected[k], atol=1e-6)
 
 
 def test_scale_voltage_current_power():
diff --git a/pvlib/tests/test_tools.py b/pvlib/tests/test_tools.py
index 42eaedca58..fa42fbf82d 100644
--- a/pvlib/tests/test_tools.py
+++ b/pvlib/tests/test_tools.py
@@ -1,6 +1,7 @@
 import pytest
 
 from pvlib import tools
+import numpy as np
 
 
 @pytest.mark.parametrize('keys, input_dict, expected', [
@@ -12,3 +13,35 @@
 def test_build_kwargs(keys, input_dict, expected):
     kwargs = tools._build_kwargs(keys, input_dict)
     assert kwargs == expected
+
+
+def _obj_test_golden_sect(params, loc):
+    return params[loc] * (1. - params['c'] * params[loc]**params['n'])
+
+
+@pytest.mark.parametrize('params, lb, ub, expected, func', [
+    ({'c': 1., 'n': 1.}, 0., 1., 0.5, _obj_test_golden_sect),
+    ({'c': 1e6, 'n': 6.}, 0., 1., 0.07230200263994839, _obj_test_golden_sect),
+    ({'c': 0.2, 'n': 0.3}, 0., 100., 89.14332727531685, _obj_test_golden_sect)
+])
+def test__golden_sect_DataFrame(params, lb, ub, expected, func):
+    v, x = tools._golden_sect_DataFrame(params, lb, ub, func)
+    assert np.isclose(x, expected, atol=1e-8)
+
+
+def test__golden_sect_DataFrame_atol():
+    params = {'c': 0.2, 'n': 0.3}
+    expected = 89.14332727531685
+    v, x = tools._golden_sect_DataFrame(
+        params, 0., 100., _obj_test_golden_sect, atol=1e-12)
+    assert np.isclose(x, expected, atol=1e-12)
+
+
+def test__golden_sect_DataFrame_vector():
+    params = {'c': np.array([1., 2.]), 'n': np.array([1., 1.])}
+    lower = np.array([0., 0.001])
+    upper = np.array([1.1, 1.2])
+    expected = np.array([0.5, 0.25])
+    v, x = tools._golden_sect_DataFrame(params, lower, upper,
+                                        _obj_test_golden_sect)
+    assert np.allclose(x, expected, atol=1e-8)
diff --git a/pvlib/tools.py b/pvlib/tools.py
index eef80a3b37..2177c75ba5 100644
--- a/pvlib/tools.py
+++ b/pvlib/tools.py
@@ -277,51 +277,61 @@ def _build_args(keys, input_dict, dict_name):
 
 # Created April,2014
 # Author: Rob Andrews, Calama Consulting
-
-def _golden_sect_DataFrame(params, VL, VH, func):
+# Modified: November, 2020 by C. W. Hansen, to add atol and change exit
+# criteria
+def _golden_sect_DataFrame(params, lower, upper, func, atol=1e-8):
     """
-    Vectorized golden section search for finding MPP from a dataframe
-    timeseries.
+    Vectorized golden section search for finding maximum of a function of a
+    single variable.
 
     Parameters
     ----------
-    params : dict
-        Dictionary containing scalars or arrays
-        of inputs to the function to be optimized.
-        Each row should represent an independent optimization.
+    params : dict or Dataframe
+        Parameters to be passed to `func`.
 
-    VL: float
-        Lower bound of the optimization
+    lower: numeric
+        Lower bound for the optimization
 
-    VH: float
-        Upper bound of the optimization
+    upper: numeric
+        Upper bound for the optimization
 
     func: function
-        Function to be optimized must be in the form f(array-like, x)
+        Function to be optimized. Must be in the form
+        result = f(dict or DataFrame, str), where result is a dict or DataFrame
+        that also contains the function output, and str is the key
+        corresponding to the function's input variable.
 
     Returns
     -------
-    func(df,'V1') : DataFrame
-        function evaluated at the optimal point
+    numeric
+        function evaluated at the optimal points
 
-    df['V1']: Dataframe
-        Dataframe of optimal points
+    numeric
+        optimal points
 
     Notes
     -----
-    This function will find the MAXIMUM of a function
+    This function will find the points where the function is maximized.
+
+    See also
+    --------
+    pvlib.singlediode._pwr_optfcn
     """
 
+    phim1 = (np.sqrt(5) - 1) / 2
+
     df = params
-    df['VH'] = VH
-    df['VL'] = VL
+    df['VH'] = upper
+    df['VL'] = lower
 
-    errflag = True
+    converged = False
     iterations = 0
+    iterlimit = 1 + np.max(
+        np.trunc(np.log(atol / (df['VH'] - df['VL'])) / np.log(phim1)))
 
-    while errflag:
+    while not converged and (iterations < iterlimit):
 
-        phi = (np.sqrt(5)-1)/2*(df['VH']-df['VL'])
+        phi = phim1 * (df['VH'] - df['VL'])
         df['V1'] = df['VL'] + phi
         df['V2'] = df['VH'] - phi
 
@@ -332,15 +342,15 @@ def _golden_sect_DataFrame(params, VL, VH, func):
         df['VL'] = df['V2']*df['SW_Flag'] + df['VL']*(~df['SW_Flag'])
         df['VH'] = df['V1']*~df['SW_Flag'] + df['VH']*(df['SW_Flag'])
 
-        err = df['V1'] - df['V2']
-        try:
-            errflag = (abs(err) > .01).any()
-        except ValueError:
-            errflag = (abs(err) > .01)
+        err = abs(df['V2'] - df['V1'])
 
+        # works with single value because err is np.float64
+        converged = (err < atol).all()
+        # err will be less than atol before iterations hit the limit
+        # but just to be safe
         iterations += 1
 
-        if iterations > 50:
-            raise Exception("EXCEPTION:iterations exceeded maximum (50)")
+    if iterations > iterlimit:
+        raise Exception("iterations exceeded maximum")  # pragma: no cover
 
     return func(df, 'V1'), df['V1']