diff --git a/README.md b/README.md
index b6b78e9..cfb9ec4 100644
--- a/README.md
+++ b/README.md
@@ -166,9 +166,14 @@ As noted above, use list comprehension if you need vector data without Numpy.
 
 ## Compare to Matlab Mapping and Aerospace Toolbox
 
-In the
-[./scripts/matlab/](scripts/matlab/)
-directory are "compare_*.py" Python scripts that use
+The tests in files tests/test_matlab*.py selected by
+
+```sh
+pytest -k matlab
+# run from pymap3d/ top-level directory
+```
+
+use
 [Matlab Engine for Python](https://www.mathworks.com/help/matlab/matlab_external/install-the-matlab-engine-for-python.html)
 to compare Python PyMap3D output with Matlab output using Matlab functions.
 
diff --git a/scripts/matlab/compare_ecef2eci.py b/scripts/matlab/compare_ecef2eci.py
deleted file mode 100644
index 8aeb856..0000000
--- a/scripts/matlab/compare_ecef2eci.py
+++ /dev/null
@@ -1,81 +0,0 @@
-#!/usr/bin/env python3
-"""
-Compare ecef2eci() with Matlab Aerospace Toolbox
-"""
-
-from __future__ import annotations
-import argparse
-import logging
-from datetime import datetime
-
-import numpy as np
-
-import pymap3d
-
-from matlab_engine import matlab_engine, has_aerospace, pydt2matdt, mat2np
-
-
-def ecef2eci(eng, utc_m, ecef):
-    if has_aerospace(eng):
-        eci = eng.ecef2eci(utc_m, np.asarray(ecef), nargout=1)
-        return mat2np(eci)
-
-    return eng.matmap3d.ecef2eci(utc_m, *ecef, nargout=3)
-
-
-def eci2ecef(eng, utc_m, eci):
-    if has_aerospace(eng):
-        ecef = eng.eci2ecef(utc_m, eci, nargout=1)
-        return mat2np(ecef)
-
-    return eng.matmap3d.eci2ecef(utc_m, *eci, nargout=3)
-
-
-def compare_ecef2eci(eng, ecef, utc: datetime) -> bool:
-    eci_py = pymap3d.ecef2eci(ecef[0], ecef[1], ecef[2], utc)
-
-    eci_m = ecef2eci(eng, pydt2matdt(eng, utc), ecef)
-
-    ok = bool(np.isclose(eci_py, eci_m, rtol=0.01).all())
-
-    if not ok:
-        logging.error(f"eci2ecef did not match Matlab  {eci_py}  !=  {eci_m}")
-
-    return ok
-
-
-def compare_eci2ecef(eng, eci, utc) -> bool:
-    ecef_py = pymap3d.eci2ecef(eci[0], eci[1], eci[2], utc)
-
-    ecef_m = eci2ecef(eng, pydt2matdt(eng, utc), eci)
-
-    ok = bool(np.isclose(ecef_py, ecef_m, rtol=0.02).all())
-
-    if not ok:
-        logging.error(f"eci2ecef did not match Matlab  {ecef_py}  !=  {ecef_m}")
-
-    return ok
-
-
-ecef = [-5762640.0, -1682738.0, 3156028.0]
-eci = [-3009680.518620539, 5194367.153184303, 3156028.0]
-utc_py = datetime(2019, 1, 4, 12)
-
-p = argparse.ArgumentParser(description="compare ecef eci")
-p.add_argument(
-    "-f", "--force_matmap3d", help="use matmap3d instead of Matlab OEM Toolbox", action="store_true"
-)
-P = p.parse_args()
-
-eng = matlab_engine()
-
-print("Mapping Toolbox:", has_aerospace(eng, P.force_matmap3d))
-
-if ecef2eci_ok := compare_ecef2eci(eng, ecef, utc_py):
-    print("OK: PyMap3d ecef2eci vs. Matlab ecef2eci")
-
-if eci2ecef_ok := compare_eci2ecef(eng, eci, utc_py):
-    print("OK: PyMap3d eci2ecef vs. Matlab eci2ecef")
-
-if not ecef2eci_ok and eci2ecef_ok:
-    raise ValueError("compare_ecef2eci: Matlab compare mismatch")
diff --git a/scripts/matlab/compare_lox.py b/scripts/matlab/compare_lox.py
deleted file mode 100644
index f14a279..0000000
--- a/scripts/matlab/compare_lox.py
+++ /dev/null
@@ -1,58 +0,0 @@
-#!/usr/bin/env python3
-"""
-Compare with Matlab Mapping toolbox reckon()
-"""
-
-from __future__ import annotations
-import argparse
-import logging
-from math import isclose
-
-from matlab_engine import matlab_engine, has_mapping
-
-from pymap3d.lox import loxodrome_direct
-
-
-def reckon(eng, lat1: float, lon1: float, rng: float, az: float) -> tuple[float, float]:
-    """Using Matlab Engine to do same thing as Pymap3d"""
-
-    if has_mapping(eng):
-        return eng.reckon("rh", lat1, lon1, rng, az, eng.wgs84Ellipsoid(), nargout=2)
-    else:
-        return eng.matmap3d.vreckon(lat1, lon1, rng, az, nargout=2)
-
-
-def stability(eng) -> bool:
-    clat, clon, rng = 35.0, 140.0, 50000.0  # arbitrary
-    ok = True
-    for i in range(20):
-        for azi in (90 + 10.0 ** (-i), -90 + 10.0 ** (-i), 270 + 10.0 ** (-i), -270 + 10.0 ** (-i)):
-            lat, lon = loxodrome_direct(clat, clon, rng, azi)
-
-            lat_matlab, lon_matlab = reckon(eng, clat, clon, rng, azi)
-            rstr = (
-                f"azimuth: {azi} lat,lon: Python: {lat}, {lon}  Matlab: {lat_matlab}, {lon_matlab}"
-            )
-            if not (
-                isclose(lat_matlab, lat, rel_tol=0.005) and isclose(lon_matlab, lon, rel_tol=0.001)
-            ):
-                ok = False
-                logging.error(rstr)
-
-    return ok
-
-
-p = argparse.ArgumentParser(description="compare reckon loxodrome")
-p.add_argument(
-    "-f", "--force_matmap3d", help="use matmap3d instead of Matlab OEM Toolbox", action="store_true"
-)
-P = p.parse_args()
-
-eng = matlab_engine()
-
-print("Mapping Toolbox:", has_mapping(eng, P.force_matmap3d))
-
-if stability(eng):
-    print("OK: lox_stability: comparison")
-else:
-    raise ValueError("FAIL: lox_stability comparison")
diff --git a/scripts/matlab/compare_vdist.py b/scripts/matlab/compare_vdist.py
deleted file mode 100644
index 00d9ceb..0000000
--- a/scripts/matlab/compare_vdist.py
+++ /dev/null
@@ -1,69 +0,0 @@
-#!/usr/bin/env python3
-"""
-Compare with Matlab Mapping Toolbox distance()
-"""
-
-from __future__ import annotations
-import argparse
-import logging
-from math import isclose, nan
-import numpy as np
-
-from matlab_engine import matlab_engine, has_mapping
-
-from pymap3d.vincenty import vdist
-
-
-def distance(eng, lat1, lon1, lat2, lon2) -> tuple[float, float]:
-    """Using Matlab Engine to do same thing as Pymap3d"""
-
-    if has_mapping(eng):
-        return eng.distance(lat1, lon1, lat2, lon2, eng.wgs84Ellipsoid(), nargout=2)
-    else:
-        return eng.matmap3d.vdist(lat1, lon1, lat2, lon2, nargout=2)
-
-
-def stability(eng) -> bool:
-    dlast, alast = nan, nan
-    lon1, lon2 = 0.0, 1.0
-    ok = True
-    for i in range(20):
-        lat1 = lat2 = 10.0 ** (-i)
-
-        dist_m, az_deg = vdist(lat1, lon1, lat2, lon2)
-
-        assert dist_m != dlast
-        assert az_deg != alast
-        dist_matlab, az_matlab = distance(eng, lat1, lon1, lat2, lon2)
-
-        assert np.isreal(dist_m), f"Python vdist distance is not real for input latitude {lat1}"
-        assert np.isreal(dist_matlab), f"Matlab distance is not real for input latitude {lat1}"
-
-        if isclose(dist_matlab, dist_m) and isclose(az_matlab, az_deg):
-            print(f"latitudes {lat1} {lat2}: {dist_m} meters {az_deg} deg azimuth")
-        else:
-            ok = False
-            logging.error(
-                f"""MISMATCH: latitude {lat1} {lat2}
-    azimuth: Python: {az_matlab}  Matlab: {az_deg}
-    distance: Python: {dist_m}  Matlab: {dist_matlab}
-    """
-            )
-
-    return ok
-
-
-p = argparse.ArgumentParser(description="compare vdist")
-p.add_argument(
-    "-f", "--force_matmap3d", help="use matmap3d instead of Matlab OEM Toolbox", action="store_true"
-)
-P = p.parse_args()
-
-eng = matlab_engine()
-
-print("Mapping Toolbox:", has_mapping(eng, P.force_matmap3d))
-
-if stability(eng):
-    print("OK: vdist compare")
-else:
-    raise ValueError("FAIL: vdist compare")
diff --git a/scripts/matlab/compare_vreckon.py b/scripts/matlab/compare_vreckon.py
deleted file mode 100644
index a557b83..0000000
--- a/scripts/matlab/compare_vreckon.py
+++ /dev/null
@@ -1,96 +0,0 @@
-#!/usr/bin/env python3
-"""
-Compare with Matlab Mapping Toolbox reckon()
-"""
-
-from __future__ import annotations
-import argparse
-import logging
-from math import isclose, nan
-import numpy as np
-
-from matlab_engine import matlab_engine, has_mapping
-
-import pymap3d.vincenty
-
-
-def reckon(eng, lat1: float, lon1: float, srng: float, az: float) -> tuple[float, float]:
-    """Using Matlab Engine to do same thing as Pymap3d"""
-
-    if has_mapping(eng):
-        return eng.reckon("gc", lat1, lon1, srng, az, eng.wgs84Ellipsoid(), nargout=2)
-    else:
-        return eng.matmap3d.vreckon(lat1, lon1, srng, az, nargout=2)
-
-
-def stability(eng) -> bool:
-    dlast, alast = nan, nan
-    lon1, lon2 = 0.0, 1.0
-    ok = True
-    for i in range(20):
-        lat1 = lat2 = 10.0 ** (-i)
-
-        dist_m, az_deg = pymap3d.vincenty.vreckon(lat1, lon1, lat2, lon2)
-
-        assert dist_m != dlast
-        assert az_deg != alast
-        dist_matlab, az_matlab = reckon(eng, lat1, lon1, lat2, lon2)
-
-        assert np.isreal(dist_m), f"Python reckon is not real for input latitude {lat1}"
-        assert np.isreal(dist_matlab), f"Matlab reckon is not real for input latitude {lat1}"
-
-        if isclose(dist_matlab, dist_m) and isclose(az_matlab, az_deg, rel_tol=0.005):
-            print(f"latitudes {lat1} {lat2}: {dist_m} meters {az_deg} deg azimuth")
-        else:
-            ok = False
-            logging.error(
-                f"""MISMATCH: latitude {lat1} {lat2}
-    azimuth: Python: {az_matlab}  Matlab: {az_deg}
-    distance: Python: {dist_m}  Matlab: {dist_matlab}
-    """
-            )
-
-    return ok
-
-
-def unit(eng) -> bool:
-    """
-    Test various extrema and other values of interest
-    """
-
-    latlon88 = 52.22610277777778, -1.2696583333333333
-    srng88 = 839.63
-    az88 = 63.02
-
-    # issue 88
-    lat_p, lon_p = pymap3d.vincenty.vreckon(*latlon88, srng88, az88)
-
-    lat_m, lon_m = reckon(eng, *latlon88, srng88, az88)
-
-    ok = isclose(lat_p, lat_m) and isclose(lon_p, lon_m)
-
-    if not ok:
-        logging.error(
-            f"""MISMATCH: lat/lon {latlon88[0]} {latlon88[1]}
-azimuth: Python: {az88}  Matlab: {az88}
-distance: Python: {srng88}  Matlab: {srng88}
-"""
-        )
-
-    return ok
-
-
-p = argparse.ArgumentParser(description="compare vreckon")
-p.add_argument(
-    "-f", "--force_matmap3d", help="use matmap3d instead of Matlab OEM Toolbox", action="store_true"
-)
-P = p.parse_args()
-
-eng = matlab_engine()
-
-print("Mapping Toolbox:", has_mapping(eng, P.force_matmap3d))
-
-if unit(eng) and stability(eng):
-    print("OK: vreckon compare")
-else:
-    raise ValueError("FAIL: vreckon compare")
diff --git a/scripts/matlab/matlab_engine.py b/scripts/matlab/matlab_engine.py
deleted file mode 100644
index 12380d2..0000000
--- a/scripts/matlab/matlab_engine.py
+++ /dev/null
@@ -1,73 +0,0 @@
-import functools
-from pathlib import Path
-from datetime import datetime
-import numpy as np
-
-import matlab.engine
-
-
-def matlab_engine():
-    cwd = Path(__file__).parent
-    eng = matlab.engine.start_matlab("-nojvm")
-    eng.addpath(eng.genpath(str(cwd)), nargout=0)
-    return eng
-
-
-def mat2np(m):
-    """
-    Matlab matrix to numpy array
-    """
-    return np.array(m.tomemoryview()).reshape(m.size, order="F")
-
-
-def pydt2matdt(eng, utc: datetime):
-    """
-    Python datetime.dateime to Matlab datetime
-    """
-    return eng.datetime(utc.year, utc.month, utc.day, utc.hour, utc.minute, utc.second)
-
-
-@functools.cache
-def has_matmap3d(eng) -> bool:
-    cwd = Path(__file__).parent
-    d = cwd.parents[2] / "matmap3d"
-    print(f"Looking in {d} for matmap3d")
-
-    if d.is_dir():
-        eng.addpath(str(d), nargout=0)
-        return True
-
-    return False
-
-
-@functools.cache
-def has_aerospace(
-    eng: matlab.engine.matlabengine.MatlabEngine, force_matmap3d: bool = False
-) -> bool:
-    if force_matmap3d:
-        if not has_matmap3d(eng):
-            raise EnvironmentError("did not find MatMap3d")
-        return False
-
-    if eng.matlab_toolbox()["aerospace"]:
-        return True
-
-    if not has_matmap3d(eng):
-        raise EnvironmentError(f"Matlab {eng.version()} does not have Aerospace Toolbox")
-
-    return False
-
-
-def has_mapping(eng: matlab.engine.matlabengine.MatlabEngine, force_matmap3d: bool = False) -> bool:
-    if force_matmap3d:
-        if not has_matmap3d(eng):
-            raise EnvironmentError("did not find MatMap3d")
-        return False
-
-    if eng.matlab_toolbox()["mapping"]:
-        return True
-
-    if not has_matmap3d(eng):
-        raise EnvironmentError(f"Matlab {eng.version()} does not have Mapping Toolbox")
-
-    return False
diff --git a/src/pymap3d/tests/matlab_engine.py b/src/pymap3d/tests/matlab_engine.py
new file mode 100644
index 0000000..f74ba47
--- /dev/null
+++ b/src/pymap3d/tests/matlab_engine.py
@@ -0,0 +1,46 @@
+import functools
+from pathlib import Path
+from datetime import datetime
+
+import matlab.engine
+
+
+@functools.cache
+def matlab_engine():
+    """
+    only cached because used by Pytest in multiple tests
+    """
+    cwd = Path(__file__).parent
+    eng = matlab.engine.start_matlab("-nojvm")
+    eng.addpath(eng.genpath(str(cwd)), nargout=0)
+    return eng
+
+
+def pydt2matdt(eng, utc: datetime):
+    """
+    Python datetime.dateime to Matlab datetime
+    """
+    return eng.datetime(utc.year, utc.month, utc.day, utc.hour, utc.minute, utc.second)
+
+
+@functools.cache
+def has_matmap3d(eng) -> bool:
+    cwd = Path(__file__).parent
+    d = cwd.parents[3] / "matmap3d"
+    print(f"Looking in {d} for matmap3d")
+
+    if d.is_dir():
+        eng.addpath(str(d), nargout=0)
+        return True
+
+    return False
+
+
+@functools.cache
+def has_aerospace(eng) -> bool:
+    return eng.matlab_toolbox()["aerospace"]
+
+
+@functools.cache
+def has_mapping(eng) -> bool:
+    return eng.matlab_toolbox()["mapping"]
diff --git a/scripts/matlab/matlab_toolbox.m b/src/pymap3d/tests/matlab_toolbox.m
similarity index 100%
rename from scripts/matlab/matlab_toolbox.m
rename to src/pymap3d/tests/matlab_toolbox.m
diff --git a/src/pymap3d/tests/test_matlab_ecef2eci.py b/src/pymap3d/tests/test_matlab_ecef2eci.py
new file mode 100644
index 0000000..425d53e
--- /dev/null
+++ b/src/pymap3d/tests/test_matlab_ecef2eci.py
@@ -0,0 +1,75 @@
+"""
+Compare ecef2eci() with Matlab Aerospace Toolbox
+"""
+
+from __future__ import annotations
+from datetime import datetime
+
+import pytest
+from pytest import approx
+
+import pymap3d
+
+try:
+    import numpy as np
+    from .matlab_engine import matlab_engine, has_aerospace, has_matmap3d, pydt2matdt
+except ImportError:
+    pytest.skip("Matlab Engine not found", allow_module_level=True)
+
+
+def ecef2eci(eng, matmap3d: bool, utc_m, ecef):
+    if matmap3d:
+        return eng.matmap3d.ecef2eci(utc_m, *ecef, nargout=3)
+
+    return np.array(eng.ecef2eci(utc_m, np.asarray(ecef), nargout=1)).squeeze()
+
+
+def eci2ecef(eng, matmap3d: bool, utc_m, eci):
+    if matmap3d:
+        return eng.matmap3d.eci2ecef(utc_m, *eci, nargout=3)
+
+    return np.array(eng.eci2ecef(utc_m, np.asarray(eci), nargout=1)).squeeze()
+
+
+@pytest.mark.parametrize("matmap3d", [False, True])
+def test_compare_ecef2eci(matmap3d):
+    eng = matlab_engine()
+
+    if matmap3d:
+        if not has_matmap3d(eng):
+            pytest.skip("Matmap3d not found")
+    else:
+        if not has_aerospace(eng):
+            pytest.skip("Aerospace Toolbox not found")
+
+    ecef = [-5762640.0, -1682738.0, 3156028.0]
+    utc = datetime(2019, 1, 4, 12)
+    rtol = 0.01
+
+    eci_py = pymap3d.ecef2eci(ecef[0], ecef[1], ecef[2], utc)
+
+    eci_m = ecef2eci(eng, matmap3d, pydt2matdt(eng, utc), ecef)
+
+    assert eci_py == approx(eci_m, rel=rtol)
+
+
+@pytest.mark.parametrize("matmap3d", [False, True])
+def test_compare_eci2ecef(matmap3d):
+    eng = matlab_engine()
+
+    if matmap3d:
+        if not has_matmap3d(eng):
+            pytest.skip("Matmap3d not found")
+    else:
+        if not has_aerospace(eng):
+            pytest.skip("Aerospace Toolbox not found")
+
+    eci = [-3009680.518620539, 5194367.153184303, 3156028.0]
+    utc = datetime(2019, 1, 4, 12)
+    rtol = 0.02
+
+    ecef_py = pymap3d.eci2ecef(eci[0], eci[1], eci[2], utc)
+
+    ecef_m = eci2ecef(eng, matmap3d, pydt2matdt(eng, utc), eci)
+
+    assert ecef_py == approx(ecef_m, rel=rtol)
diff --git a/src/pymap3d/tests/test_matlab_lox.py b/src/pymap3d/tests/test_matlab_lox.py
new file mode 100644
index 0000000..ad0e8d3
--- /dev/null
+++ b/src/pymap3d/tests/test_matlab_lox.py
@@ -0,0 +1,51 @@
+#!/usr/bin/env python3
+"""
+Compare with Matlab Mapping toolbox reckon()
+"""
+
+from __future__ import annotations
+
+import pytest
+from pytest import approx
+
+try:
+    from .matlab_engine import matlab_engine, has_matmap3d, has_mapping
+except ImportError:
+    pytest.skip("Matlab Engine not found", allow_module_level=True)
+
+
+from pymap3d.lox import loxodrome_direct
+
+
+def reckon(
+    eng, matmap3d: bool, lat1: float, lon1: float, rng: float, az: float
+) -> tuple[float, float]:
+    """Using Matlab Engine to do same thing as Pymap3d"""
+
+    if matmap3d:
+        return eng.matmap3d.vreckon(lat1, lon1, rng, az, nargout=2)
+
+    return eng.reckon("rh", lat1, lon1, rng, az, eng.wgs84Ellipsoid(), nargout=2)
+
+
+@pytest.mark.parametrize("matmap3d", [False, True])
+def test_lox_stability(matmap3d):
+    eng = matlab_engine()
+
+    if matmap3d:
+        if not has_matmap3d(eng):
+            pytest.skip("Matmap3d not found")
+    else:
+        if not has_mapping(eng):
+            pytest.skip("Matlab Toolbox not found")
+
+    clat, clon, rng = 35.0, 140.0, 50000.0  # arbitrary
+
+    for i in range(20):
+        for azi in (90 + 10.0 ** (-i), -90 + 10.0 ** (-i), 270 + 10.0 ** (-i), -270 + 10.0 ** (-i)):
+            lat, lon = loxodrome_direct(clat, clon, rng, azi)
+
+            lat_matlab, lon_matlab = reckon(eng, matmap3d, clat, clon, rng, azi)
+
+            assert lat == approx(lat_matlab, rel=0.005)
+            assert lon == approx(lon_matlab, rel=0.001)
diff --git a/src/pymap3d/tests/test_matlab_vdist.py b/src/pymap3d/tests/test_matlab_vdist.py
new file mode 100644
index 0000000..c75c713
--- /dev/null
+++ b/src/pymap3d/tests/test_matlab_vdist.py
@@ -0,0 +1,52 @@
+"""
+Compare with Matlab Mapping Toolbox distance()
+"""
+
+from __future__ import annotations
+from math import nan
+
+import pytest
+from pytest import approx
+
+try:
+    from .matlab_engine import matlab_engine, has_mapping, has_matmap3d
+except ImportError:
+    pytest.skip("Matlab Engine not found", allow_module_level=True)
+
+from pymap3d.vincenty import vdist
+
+
+def distance(eng, matmap3d: bool, lat1, lon1, lat2, lon2) -> tuple[float, float]:
+    """Using Matlab Engine to do same thing as Pymap3d"""
+
+    if matmap3d:
+        return eng.matmap3d.vdist(lat1, lon1, lat2, lon2, nargout=2)
+
+    return eng.distance(lat1, lon1, lat2, lon2, eng.wgs84Ellipsoid(), nargout=2)
+
+
+@pytest.mark.parametrize("matmap3d", [False, True])
+def test_matlab_stability(matmap3d):
+    eng = matlab_engine()
+
+    if matmap3d:
+        if not has_matmap3d(eng):
+            pytest.skip("Matmap3d not found")
+    else:
+        if not has_mapping(eng):
+            pytest.skip("Matlab Toolbox not found")
+
+    dlast, alast = nan, nan
+    lon1, lon2 = 0.0, 1.0
+
+    for i in range(20):
+        lat1 = lat2 = 10.0 ** (-i)
+
+        dist_m, az_deg = vdist(lat1, lon1, lat2, lon2)
+
+        assert dist_m != dlast
+        assert az_deg != alast
+        dist_matlab, az_matlab = distance(eng, matmap3d, lat1, lon1, lat2, lon2)
+
+        assert dist_m == approx(dist_matlab)
+        assert az_deg == approx(az_matlab)
diff --git a/src/pymap3d/tests/test_matlab_vreckon.py b/src/pymap3d/tests/test_matlab_vreckon.py
new file mode 100644
index 0000000..4dbbde8
--- /dev/null
+++ b/src/pymap3d/tests/test_matlab_vreckon.py
@@ -0,0 +1,82 @@
+"""
+Compare with Matlab Mapping Toolbox reckon()
+"""
+
+from __future__ import annotations
+from math import nan
+
+import pytest
+from pytest import approx
+
+try:
+    from .matlab_engine import matlab_engine, has_mapping, has_matmap3d
+except ImportError:
+    pytest.skip("Matlab Engine not found", allow_module_level=True)
+
+import pymap3d.vincenty
+
+
+def reckon(
+    eng, matmap3d: bool, lat1: float, lon1: float, srng: float, az: float
+) -> tuple[float, float]:
+    """Using Matlab Engine to do same thing as Pymap3d"""
+
+    if matmap3d:
+        return eng.matmap3d.vreckon(lat1, lon1, srng, az, nargout=2)
+
+    return eng.reckon("gc", lat1, lon1, srng, az, eng.wgs84Ellipsoid(), nargout=2)
+
+
+@pytest.mark.parametrize("matmap3d", [False, True])
+def test_reckon_stability(matmap3d):
+    eng = matlab_engine()
+
+    if matmap3d:
+        if not has_matmap3d(eng):
+            pytest.skip("Matmap3d not found")
+    else:
+        if not has_mapping(eng):
+            pytest.skip("Matlab Toolbox not found")
+
+    dlast, alast = nan, nan
+    lon1, lon2 = 0.0, 1.0
+    for i in range(20):
+        lat1 = lat2 = 10.0 ** (-i)
+
+        dist_m, az_deg = pymap3d.vincenty.vreckon(lat1, lon1, lat2, lon2)
+
+        assert dist_m != dlast
+        assert az_deg != alast
+        dist_matlab, az_matlab = reckon(eng, matmap3d, lat1, lon1, lat2, lon2)
+
+        assert dist_m == approx(dist_matlab)
+
+        assert az_deg == approx(az_matlab, rel=0.005)
+
+
+@pytest.mark.parametrize("matmap3d", [False, True])
+def test_reckon_unit(matmap3d):
+    """
+    Test various extrema and other values of interest
+    """
+
+    eng = matlab_engine()
+
+    if matmap3d:
+        if not has_matmap3d(eng):
+            pytest.skip("Matmap3d not found")
+    else:
+        if not has_mapping(eng):
+            pytest.skip("Matlab Toolbox not found")
+
+    latlon88 = 52.22610277777778, -1.2696583333333333
+    srng88 = 839.63
+    az88 = 63.02
+
+    # issue 88
+    lat_p, lon_p = pymap3d.vincenty.vreckon(*latlon88, srng88, az88)
+
+    lat_m, lon_m = reckon(eng, matmap3d, *latlon88, srng88, az88)
+
+    assert lat_p == approx(lat_m)
+    assert lon_p == approx(lon_m)