diff --git a/docs/source/Support/bskReleaseNotes.rst b/docs/source/Support/bskReleaseNotes.rst
index ba2ca8b386..11dba24731 100644
--- a/docs/source/Support/bskReleaseNotes.rst
+++ b/docs/source/Support/bskReleaseNotes.rst
@@ -109,6 +109,7 @@ Version  |release|
 - Updated the CI build that includes the documentation to fail if a doxygen warning happens
 - Removed deprecated swig code that allowed still importing `sys_model.h` instead of `sys_model.i`
 - Updated :ref:`groundMapping` to correct behavior if ``maximumRange == -1``
+- Updated scripts to work with ``matplotlib`` version 3.10.x without errors or warnings
 
 
 Version 2.5.0 (Sept. 30, 2024)
diff --git a/examples/BskSim/plotting/BSK_Plotting.py b/examples/BskSim/plotting/BSK_Plotting.py
index 7db30405b6..d66f9f70ae 100644
--- a/examples/BskSim/plotting/BSK_Plotting.py
+++ b/examples/BskSim/plotting/BSK_Plotting.py
@@ -271,7 +271,7 @@ def plot_rw_friction(timeData, dataFrictionRW, numRW, dataFaultLog=[],  id=None,
 
 def plot_planet(oe, planet):
     b = oe.a * np.sqrt(1 - oe.e * oe.e)
-    plt.figure(figsize=np.array((1.0, b / oe.a)) * 4.75, dpi=100)
+    plt.figure(figsize=tuple(np.array((1.0, b / oe.a)) * 4.75), dpi=100)
     plt.axis(np.array([-oe.a, oe.a, -b, b]) / 1000 * 1.75)
     # draw the planet
     fig = plt.gcf()
diff --git a/examples/scenarioAttGuideHyperbolic.py b/examples/scenarioAttGuideHyperbolic.py
index 0b2b3668ed..a95a6f1524 100755
--- a/examples/scenarioAttGuideHyperbolic.py
+++ b/examples/scenarioAttGuideHyperbolic.py
@@ -173,7 +173,7 @@ def plot_orbit(oe, mu, planet_radius, dataPos, dataVel):
     """Plot the spacecraft orbit trajectory."""
     # draw orbit in perifocal frame
     p = oe.a * (1 - oe.e * oe.e)
-    plt.figure(4, figsize=np.array((1.0, 1.)) * 4.75, dpi=100)
+    plt.figure(4, figsize=tuple(np.array((1.0, 1.)) * 4.75), dpi=100)
     # draw the planet
     fig = plt.gcf()
     ax = fig.gca()
diff --git a/examples/scenarioBasicOrbit.py b/examples/scenarioBasicOrbit.py
index 78af5b50fe..c90540e66a 100644
--- a/examples/scenarioBasicOrbit.py
+++ b/examples/scenarioBasicOrbit.py
@@ -487,7 +487,7 @@ def plotOrbits(timeAxis, posData, velData, oe, mu, P, orbitCase, useSphericalHar
         # draw orbit in perifocal frame
         b = oe.a * np.sqrt(1 - oe.e * oe.e)
         p = oe.a * (1 - oe.e * oe.e)
-        plt.figure(2, figsize=np.array((1.0, b / oe.a)) * 4.75, dpi=100)
+        plt.figure(2, figsize=tuple(np.array((1.0, b / oe.a)) * 4.75), dpi=100)
         plt.axis(np.array([-oe.rApoap, oe.rPeriap, -b, b]) / 1000 * 1.25)
         # draw the planet
         fig = plt.gcf()
diff --git a/examples/scenarioBasicOrbitStream.py b/examples/scenarioBasicOrbitStream.py
index 544e75ffe6..14d7982067 100644
--- a/examples/scenarioBasicOrbitStream.py
+++ b/examples/scenarioBasicOrbitStream.py
@@ -426,7 +426,7 @@ def run(show_plots, liveStream, broadcastStream, timeStep, orbitCase, useSpheric
         # draw orbit in perifocal frame
         b = oe.a * np.sqrt(1 - oe.e * oe.e)
         p = oe.a * (1 - oe.e * oe.e)
-        plt.figure(2, figsize=np.array((1.0, b / oe.a)) * 4.75, dpi=100)
+        plt.figure(2, figsize=tuple(np.array((1.0, b / oe.a)) * 4.75), dpi=100)
         plt.axis(np.array([-oe.rApoap, oe.rPeriap, -b, b]) / 1000 * 1.25)
         # draw the planet
         fig = plt.gcf()
diff --git a/examples/scenarioCentralBody.py b/examples/scenarioCentralBody.py
index f82cc161bb..23346b44c8 100644
--- a/examples/scenarioCentralBody.py
+++ b/examples/scenarioCentralBody.py
@@ -285,7 +285,7 @@ def run(show_plots, useCentral):
     # draw orbit in perifocal frame
     b = oe.a * np.sqrt(1 - oe.e * oe.e)
     p = oe.a * (1 - oe.e * oe.e)
-    plt.figure(2, figsize=np.array((1.0, b / oe.a)) * 4.75, dpi=100)
+    plt.figure(2, figsize=tuple(np.array((1.0, b / oe.a)) * 4.75), dpi=100)
     plt.axis(np.array([-oe.rApoap, oe.rPeriap, -b, b]) / 1000 * 1.25)
     # draw the planet
     fig = plt.gcf()
diff --git a/examples/scenarioDragDeorbit.py b/examples/scenarioDragDeorbit.py
index 64d205d8d6..d6c876e5d5 100644
--- a/examples/scenarioDragDeorbit.py
+++ b/examples/scenarioDragDeorbit.py
@@ -312,7 +312,7 @@ def register_fig(i):
 
     # draw orbit in perifocal frame
     b = oe.a * np.sqrt(1 - oe.e * oe.e)
-    plt.figure(1, figsize=np.array((1.0, b / oe.a)) * 4.75, dpi=100)
+    plt.figure(1, figsize=tuple(np.array((1.0, b / oe.a)) * 4.75), dpi=100)
     plt.axis(np.array([-oe.rApoap, oe.rPeriap, -b, b]) / 1000 * 1.25)
     # draw the planet
     fig, ax = register_fig(1)
diff --git a/examples/scenarioHaloOrbit.py b/examples/scenarioHaloOrbit.py
index 52f520d40d..8efeb08802 100644
--- a/examples/scenarioHaloOrbit.py
+++ b/examples/scenarioHaloOrbit.py
@@ -204,7 +204,7 @@ def run(showPlots=True):
     b = oe.a * np.sqrt(1 - oe.e * oe.e)
 
     # First plot: Draw orbit in inertial frame
-    fig = plt.figure(1, figsize=np.array((1.0, b / oe.a)) * 4.75, dpi=100)
+    fig = plt.figure(1, figsize=tuple(np.array((1.0, b / oe.a)) * 4.75), dpi=100)
     plt.axis(np.array([-oe.rApoap, oe.rPeriap, -b, b]) / 1000 * 1.25)
     ax = fig.gca()
     ax.ticklabel_format(style='scientific', scilimits=[-5, 3])
@@ -241,7 +241,7 @@ def run(showPlots=True):
 
     # Second plot: Draw orbit in frame rotating with the Moon (the center is L2 point)
     # x axis is moon position vector direction and y axis is moon velocity vector direction
-    fig = plt.figure(2, figsize=np.array((1.0, b / oe.a)) * 4.75, dpi=100)
+    fig = plt.figure(2, figsize=tuple(np.array((1.0, b / oe.a)) * 4.75), dpi=100)
     plt.axis(np.array([-1e5, 5e5, -3e5, 3e5])  * 1.25)
     ax = fig.gca()
     ax.ticklabel_format(style='scientific', scilimits=[-5, 3])
@@ -279,7 +279,7 @@ def run(showPlots=True):
     # Third plot: Draw orbit in frame rotating with the Moon (the center is L2 point)
     # x axis is moon position vector direction and y axis is the cross product direction of the moon position vector and
     # velocity vector
-    fig = plt.figure(3, figsize=np.array((1.0, b / oe.a)) * 4.75, dpi=100)
+    fig = plt.figure(3, figsize=tuple(np.array((1.0, b / oe.a)) * 4.75), dpi=100)
     plt.axis(np.array([-1e5, 5e5, -3e5, 3e5]) * 1.25)
     ax = fig.gca()
     ax.ticklabel_format(style='scientific', scilimits=[-5, 3])
@@ -316,5 +316,3 @@ def run(showPlots=True):
     run(
         True    # Show plots
     )
-
-
diff --git a/examples/scenarioIntegrators.py b/examples/scenarioIntegrators.py
index 9380a17f85..afe67ac44b 100755
--- a/examples/scenarioIntegrators.py
+++ b/examples/scenarioIntegrators.py
@@ -328,7 +328,8 @@ def run(show_plots, integratorCase):
     # draw orbit in perifocal frame
     b = oe.a * np.sqrt(1 - oe.e * oe.e)
     p = oe.a * (1 - oe.e * oe.e)
-    plt.figure(1, figsize=np.array((1.0, b / oe.a)) * 4.75, dpi=100)
+    if integratorCase == 'rk4':
+        plt.figure(1, figsize=tuple(np.array((1.0, b / oe.a)) * 4.75), dpi=100)
     plt.axis(np.array([-oe.rApoap, oe.rPeriap, -b, b]) / 1000 * 1.25)
     # draw the planet
     fig = plt.gcf()
diff --git a/examples/scenarioLagrangePointOrbit.py b/examples/scenarioLagrangePointOrbit.py
index 569b3bae02..3797f1ba7f 100644
--- a/examples/scenarioLagrangePointOrbit.py
+++ b/examples/scenarioLagrangePointOrbit.py
@@ -273,7 +273,7 @@ def run(lagrangePoint, nOrbits, timestep, showPlots=True):
     b = oe.a * np.sqrt(1 - oe.e * oe.e)
 
     # First plot: Draw orbit in inertial frame
-    fig = plt.figure(1, figsize=np.array((1.0, b / oe.a)) * 4.75, dpi=100)
+    fig = plt.figure(1, figsize=tuple(np.array((1.0, b / oe.a)) * 4.75), dpi=100)
     plt.axis(np.array([-oe.rApoap, oe.rPeriap, -b, b]) / 1000 * 1.25)
     ax = fig.gca()
     ax.ticklabel_format(style='scientific', scilimits=[-5, 3])
@@ -317,7 +317,7 @@ def run(lagrangePoint, nOrbits, timestep, showPlots=True):
     figureList[pltName] = plt.figure(1)
 
     # Second plot: Draw orbit in frame rotating with the Moon
-    fig = plt.figure(2, figsize=np.array((1.0, b / oe.a)) * 4.75, dpi=100)
+    fig = plt.figure(2, figsize=tuple(np.array((1.0, b / oe.a)) * 4.75), dpi=100)
     plt.axis(np.array([-oe.rApoap, oe.rPeriap, -b, b]) / 1000 * 1.25)
     ax = fig.gca()
     ax.ticklabel_format(style='scientific', scilimits=[-5, 3])
@@ -387,5 +387,3 @@ def run(lagrangePoint, nOrbits, timestep, showPlots=True):
         300,       # Timestep (seconds)
         True    # Show plots
     )
-
-
diff --git a/examples/scenarioOrbitMultiBody.py b/examples/scenarioOrbitMultiBody.py
index 7248d9d16f..9a7a8ed7b9 100755
--- a/examples/scenarioOrbitMultiBody.py
+++ b/examples/scenarioOrbitMultiBody.py
@@ -336,7 +336,7 @@ def run(show_plots, scCase):
         omega0 = oeData.omega
         b = oeData.a * np.sqrt(1 - oeData.e * oeData.e)
         p = oeData.a * (1 - oeData.e * oeData.e)
-        plt.figure(2, figsize=np.array((1.0, b / oeData.a)) * 4.75, dpi=100)
+        plt.figure(2, figsize=tuple(np.array((1.0, b / oeData.a)) * 4.75), dpi=100)
         plt.axis(np.array([-oeData.rApoap, oeData.rPeriap, -b, b]) / 1000 * 1.25)
 
         # draw the planet
diff --git a/examples/scenarioPatchedConics.py b/examples/scenarioPatchedConics.py
index ed555eee04..47d68c72e0 100644
--- a/examples/scenarioPatchedConics.py
+++ b/examples/scenarioPatchedConics.py
@@ -393,7 +393,7 @@ def run(show_plots):
     plt.close("all")  # clears out plots from earlier test runs
     b = oe.a * np.sqrt(1 - oe.e * oe.e)
     p = oe.a * (1 - oe.e * oe.e)
-    plt.figure(1) #, figsize=np.array((1.0, b / oe.a)) * 4.75, dpi=100)
+    plt.figure(1) #, figsize=tuple(np.array((1.0, b / oe.a)) * 4.75), dpi=100)
     # plt.axis(np.array([-oe.rApoap, oe.rPeriap, -b, b]) / 1000 * 1.25)
     plt.axis('equal')
     plt.axis([-20000, 50000, -10000, 10000])
diff --git a/requirements.txt b/requirements.txt
index 05436dd848..4613b9fd10 100644
--- a/requirements.txt
+++ b/requirements.txt
@@ -1,5 +1,5 @@
 pandas
-matplotlib<3.10.0
+matplotlib
 numpy
 colorama
 tqdm
diff --git a/src/fswAlgorithms/smallBodyNavigation/smallBodyNavEKF/_UnitTest/test_smallBodyNavEKF.py b/src/fswAlgorithms/smallBodyNavigation/smallBodyNavEKF/_UnitTest/test_smallBodyNavEKF.py
index e17ba7f9d6..7f6a602c16 100644
--- a/src/fswAlgorithms/smallBodyNavigation/smallBodyNavEKF/_UnitTest/test_smallBodyNavEKF.py
+++ b/src/fswAlgorithms/smallBodyNavigation/smallBodyNavEKF/_UnitTest/test_smallBodyNavEKF.py
@@ -1,12 +1,12 @@
-# 
+#
 #  ISC License
-# 
+#
 #  Copyright (c) 2021, Autonomous Vehicle Systems Lab, University of Colorado Boulder
-# 
+#
 #  Permission to use, copy, modify, and/or distribute this software for any
 #  purpose with or without fee is hereby granted, provided that the above
 #  copyright notice and this permission notice appear in all copies.
-# 
+#
 #  THE SOFTWARE IS PROVIDED "AS IS" AND THE AUTHOR DISCLAIMS ALL WARRANTIES
 #  WITH REGARD TO THIS SOFTWARE INCLUDING ALL IMPLIED WARRANTIES OF
 #  MERCHANTABILITY AND FITNESS. IN NO EVENT SHALL THE AUTHOR BE LIABLE FOR
@@ -14,8 +14,8 @@
 #  WHATSOEVER RESULTING FROM LOSS OF USE, DATA OR PROFITS, WHETHER IN AN
 #  ACTION OF CONTRACT, NEGLIGENCE OR OTHER TORTIOUS ACTION, ARISING OUT OF
 #  OR IN CONNECTION WITH THE USE OR PERFORMANCE OF THIS SOFTWARE.
-# 
-# 
+#
+#
 
 import numpy as np
 from Basilisk.architecture import messaging
@@ -136,8 +136,8 @@ def smallBodyNavEKFTestFunction(show_plots):
         true_x_hat, np.array([x_hat_c_wrapped[-1,:]]), 0.1, "x_hat_c_wrapped",
         testFailCount, testMessages)
 
-    plt.figure(1)
-    plt.clf()
+    plt.close('all')
+
     plt.figure(1, figsize=(7, 5), dpi=80, facecolor='w', edgecolor='k')
     plt.ticklabel_format(useOffset=False)
     plt.plot(navTransOutMsgRec.times() * 1.0E-9, x_hat[:,0], label='x-pos')
@@ -148,8 +148,6 @@ def smallBodyNavEKFTestFunction(show_plots):
     plt.ylabel('r_BO_O (m)')
     plt.title('Estimated Relative Spacecraft Position')
 
-    plt.figure(2)
-    plt.clf()
     plt.figure(2, figsize=(7, 5), dpi=80, facecolor='w', edgecolor='k')
     plt.plot(navTransOutMsgRec.times() * 1.0E-9, x_hat[:,3], label='x-vel')
     plt.plot(navTransOutMsgRec.times() * 1.0E-9, x_hat[:,4], label='y-vel')
@@ -159,8 +157,6 @@ def smallBodyNavEKFTestFunction(show_plots):
     plt.ylabel('v_BO_O (m/s)')
     plt.title('Estimated Spacecraft Velocity')
 
-    plt.figure(5)
-    plt.clf()
     plt.figure(5, figsize=(7, 5), dpi=80, facecolor='w', edgecolor='k')
     plt.plot(navTransOutMsgRec.times() * 1.0E-9, x_hat[:,6], label='s1')
     plt.plot(navTransOutMsgRec.times() * 1.0E-9, x_hat[:,7], label='s2')
@@ -170,8 +166,6 @@ def smallBodyNavEKFTestFunction(show_plots):
     plt.ylabel('sigma_AN (rad)')
     plt.title('Estimated Asteroid Attitude')
 
-    plt.figure(6)
-    plt.clf()
     plt.figure(6, figsize=(7, 5), dpi=80, facecolor='w', edgecolor='k')
     plt.plot(navTransOutMsgRec.times() * 1.0E-9, x_hat[:,9], label='omega1')
     plt.plot(navTransOutMsgRec.times() * 1.0E-9, x_hat[:,10], label='omega2')
diff --git a/src/fswAlgorithms/smallBodyNavigation/smallBodyNavUKF/_UnitTest/test_smallBodyNavUKF.py b/src/fswAlgorithms/smallBodyNavigation/smallBodyNavUKF/_UnitTest/test_smallBodyNavUKF.py
index 009ba44fa6..1394aa4708 100644
--- a/src/fswAlgorithms/smallBodyNavigation/smallBodyNavUKF/_UnitTest/test_smallBodyNavUKF.py
+++ b/src/fswAlgorithms/smallBodyNavigation/smallBodyNavUKF/_UnitTest/test_smallBodyNavUKF.py
@@ -1,12 +1,12 @@
-# 
+#
 #  ISC License
-# 
+#
 #  Copyright (c) 2021, Autonomous Vehicle Systems Lab, University of Colorado Boulder
-# 
+#
 #  Permission to use, copy, modify, and/or distribute this software for any
 #  purpose with or without fee is hereby granted, provided that the above
 #  copyright notice and this permission notice appear in all copies.
-# 
+#
 #  THE SOFTWARE IS PROVIDED "AS IS" AND THE AUTHOR DISCLAIMS ALL WARRANTIES
 #  WITH REGARD TO THIS SOFTWARE INCLUDING ALL IMPLIED WARRANTIES OF
 #  MERCHANTABILITY AND FITNESS. IN NO EVENT SHALL THE AUTHOR BE LIABLE FOR
@@ -14,8 +14,8 @@
 #  WHATSOEVER RESULTING FROM LOSS OF USE, DATA OR PROFITS, WHETHER IN AN
 #  ACTION OF CONTRACT, NEGLIGENCE OR OTHER TORTIOUS ACTION, ARISING OUT OF
 #  OR IN CONNECTION WITH THE USE OR PERFORMANCE OF THIS SOFTWARE.
-# 
-# 
+#
+#
 
 import numpy as np
 from Basilisk.architecture import messaging
@@ -134,8 +134,7 @@ def smallBodyNavUKFTestFunction(show_plots):
         [true_x_hat], np.array([x_hat_c_wrapped[-1,:]]), 0.01, "x_hat_c_wrapped",
         testFailCount, testMessages)
 
-    plt.figure(1)
-    plt.clf()
+    plt.close('all')
     plt.figure(1, figsize=(7, 5), dpi=80, facecolor='w', edgecolor='k')
     plt.ticklabel_format(useOffset=False)
     plt.plot(smallBodyNavUKFOutMsgRec.times() * 1.0E-9 / 60, x_hat[:,0] / 1000, label='x-pos')
@@ -146,8 +145,6 @@ def smallBodyNavUKFTestFunction(show_plots):
     plt.ylabel('${}^{A}r_{BA}$ (km)')
     plt.title('Estimated Relative Spacecraft Position')
 
-    plt.figure(2)
-    plt.clf()
     plt.figure(2, figsize=(7, 5), dpi=80, facecolor='w', edgecolor='k')
     plt.plot(smallBodyNavUKFOutMsgRec.times() * 1.0E-9 / 60, x_hat[:,3], label='x-vel')
     plt.plot(smallBodyNavUKFOutMsgRec.times() * 1.0E-9 / 60, x_hat[:,4], label='y-vel')
@@ -157,8 +154,6 @@ def smallBodyNavUKFTestFunction(show_plots):
     plt.ylabel('${}^{A}v_{BA}$ (m/s)')
     plt.title('Estimated Spacecraft Velocity')
 
-    plt.figure(3)
-    plt.clf()
     plt.figure(3, figsize=(7, 5), dpi=80, facecolor='w', edgecolor='k')
     plt.plot(smallBodyNavUKFOutMsgRec.times() * 1.0E-9 / 60, x_hat[:,6], label='x-acc')
     plt.plot(smallBodyNavUKFOutMsgRec.times() * 1.0E-9 / 60, x_hat[:,7], label='y-acc')
diff --git a/src/simulation/dynamics/DynOutput/orbElemConvert/_UnitTest/test_orb_elem_convert.py b/src/simulation/dynamics/DynOutput/orbElemConvert/_UnitTest/test_orb_elem_convert.py
index f5ec847b75..dc5d90768b 100644
--- a/src/simulation/dynamics/DynOutput/orbElemConvert/_UnitTest/test_orb_elem_convert.py
+++ b/src/simulation/dynamics/DynOutput/orbElemConvert/_UnitTest/test_orb_elem_convert.py
@@ -474,9 +474,9 @@ def orbElem(a, e, i, AN, AP, f, mu, name, DispPlot):
     # create plot
     # txt = 'e = ' + str(e) + ' and a = ' + str(a) + 'km'
     fact = (len(str(abs(a)))-3.0)
+    plt.close('all')
 
     plt.figure(1,figsize=(7, 5), dpi=80, facecolor='w', edgecolor='k')
-    plt.clf()
     # fig1.text(.5, .05, txt, ha='center')
     ax1 = plt.subplot(211)
     ax1.cla()
diff --git a/src/simulation/dynamics/Integrators/_UnitTest/test_Integrators.py b/src/simulation/dynamics/Integrators/_UnitTest/test_Integrators.py
index 543fc33459..ce3aba0715 100755
--- a/src/simulation/dynamics/Integrators/_UnitTest/test_Integrators.py
+++ b/src/simulation/dynamics/Integrators/_UnitTest/test_Integrators.py
@@ -204,13 +204,14 @@ def run(doUnitTests, show_plots, integratorCase):
     #
     np.set_printoptions(precision=16)
     fileNameString = filename[len(path) + 6:-3]
-    if integratorCase == "rk4":
+    if integratorCase == "bogackiShampine":
         plt.close("all")  # clears out plots from earlier test runs
 
     # draw orbit in perifocal frame
     b = oe.a * np.sqrt(1 - oe.e * oe.e)
     p = oe.a * (1 - oe.e * oe.e)
-    plt.figure(1, figsize=np.array((1.0, b / oe.a)) * 4.75, dpi=100)
+    if integratorCase == "rk4":
+        plt.figure(1, figsize=tuple(np.array((1.0, b / oe.a)) * 4.75), dpi=100)
     plt.axis(np.array([-oe.rApoap, oe.rPeriap, -b, b]) / 1000 * 1.25)
     # draw the planet
     fig = plt.gcf()
diff --git a/src/simulation/dynamics/dragEffector/_UnitTest/test_atmoDrag.py b/src/simulation/dynamics/dragEffector/_UnitTest/test_atmoDrag.py
index 49f27c0330..038383512d 100644
--- a/src/simulation/dynamics/dragEffector/_UnitTest/test_atmoDrag.py
+++ b/src/simulation/dynamics/dragEffector/_UnitTest/test_atmoDrag.py
@@ -301,7 +301,7 @@ def run(show_plots, orbitCase, planetCase):
         # draw orbit in perifocal frame
         b = oe.a*np.sqrt(1-oe.e*oe.e)
         p = oe.a*(1-oe.e*oe.e)
-        plt.figure(2,figsize=np.array((1.0, b/oe.a))*4.75,dpi=100)
+        plt.figure(2,figsize=tuple(np.array((1.0, b/oe.a))*4.75),dpi=100)
         plt.axis(np.array([-oe.rApoap, oe.rPeriap, -b, b])/1000*1.25)
         # draw the planet
         fig = plt.gcf()
diff --git a/src/simulation/navigation/planetNav/_UnitTest/test_planetNav.py b/src/simulation/navigation/planetNav/_UnitTest/test_planetNav.py
index b637e87b91..3a165f6818 100644
--- a/src/simulation/navigation/planetNav/_UnitTest/test_planetNav.py
+++ b/src/simulation/navigation/planetNav/_UnitTest/test_planetNav.py
@@ -198,9 +198,7 @@ def planetNavTestFunction(show_plots):
         if count < 1:
             testFailCount += 1
             testMessages.append("FAILED: Too few error counts - " + str(count))
-
-    plt.figure(1)
-    plt.clf()
+    plt.close('all')
     plt.figure(1, figsize=(7, 5), dpi=80, facecolor='w', edgecolor='k')
     plt.plot(ephemerisOutMsgRec.times() * 1.0E-9, r_BN_N[:,0], label='x-position')
     plt.plot(ephemerisOutMsgRec.times() * 1.0E-9, r_BN_N[:,1], label='y-position')
@@ -210,8 +208,6 @@ def planetNavTestFunction(show_plots):
     plt.xlabel('Time (s)')
     plt.ylabel('Position (m)')
 
-    plt.figure(2)
-    plt.clf()
     plt.figure(2, figsize=(7, 5), dpi=80, facecolor='w', edgecolor='k')
     plt.plot(ephemerisOutMsgRec.times() * 1.0E-9, v_BN_N[:,0], label='x-velocity')
     plt.plot(ephemerisOutMsgRec.times() * 1.0E-9, v_BN_N[:,1], label='y-velocity')
@@ -221,8 +217,7 @@ def planetNavTestFunction(show_plots):
     plt.xlabel('Time (s)')
     plt.ylabel('Velocity (m/s)')
 
-    plt.figure(3)
-    plt.clf()
+
     plt.figure(3, figsize=(7, 5), dpi=80, facecolor='w', edgecolor='k')
     plt.plot(ephemerisOutMsgRec.times() * 1.0E-9, sigma_BN[:, 0], label='x-rotation')
     plt.plot(ephemerisOutMsgRec.times() * 1.0E-9, sigma_BN[:, 1], label='y-rotation')
@@ -232,8 +227,6 @@ def planetNavTestFunction(show_plots):
     plt.xlabel('Time (s)')
     plt.ylabel('Attitude (rad)')
 
-    plt.figure(4)
-    plt.clf()
     plt.figure(4, figsize=(7, 5), dpi=80, facecolor='w', edgecolor='k')
     plt.plot(ephemerisOutMsgRec.times() * 1.0E-9, omega_BN_B[:, 0], label='x-angular vel.')
     plt.plot(ephemerisOutMsgRec.times() * 1.0E-9, omega_BN_B[:, 1], label='y-angular vel.')
@@ -245,7 +238,7 @@ def planetNavTestFunction(show_plots):
 
     if show_plots:
         plt.show()
-        plt.close('all')
+    plt.close('all')
 
     # Corner case usage
     pMatrixBad = [[0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.],
diff --git a/src/simulation/navigation/simpleNav/_UnitTest/test_simpleNav.py b/src/simulation/navigation/simpleNav/_UnitTest/test_simpleNav.py
index d9017a3901..ac3c5d3708 100755
--- a/src/simulation/navigation/simpleNav/_UnitTest/test_simpleNav.py
+++ b/src/simulation/navigation/simpleNav/_UnitTest/test_simpleNav.py
@@ -232,8 +232,6 @@ def unitSimpleNav(show_plots):
             testFailCount += 1
             testMessages.append("FAILED: Too few error counts -" + str(count))
 
-    plt.figure(1)
-    plt.clf()
     plt.figure(1, figsize=(7, 5), dpi=80, facecolor='w', edgecolor='k')
     plt.plot(dataTransLog.times() * 1.0E-9, posNav[:,0], label='x-position')
     plt.plot(dataTransLog.times() * 1.0E-9, posNav[:,1], label='y-position')
@@ -247,8 +245,6 @@ def unitSimpleNav(show_plots):
         plt.show()
         plt.close('all')
 
-    plt.figure(2)
-    plt.clf()
     plt.figure(2, figsize=(7, 5), dpi=80, facecolor='w', edgecolor='k')
     plt.plot(dataAttLog.times() * 1.0E-9, attNav[:, 0], label='x-rotation')
     plt.plot(dataAttLog.times() * 1.0E-9, attNav[:, 1], label='y-rotation')
@@ -260,7 +256,7 @@ def unitSimpleNav(show_plots):
     unitTestSupport.writeFigureLaTeX('SimpleNavAtt', 'Simple Navigation Att Signal', plt, r'height=0.4\textwidth, keepaspectratio', path)
     if show_plots:
         plt.show()
-        plt.close('all')
+    plt.close('all')
 
     # Corner case usage
     pMatrixBad = [[0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.],
diff --git a/src/simulation/sensors/coarseSunSensor/_UnitTest/test_coarseSunSensor.py b/src/simulation/sensors/coarseSunSensor/_UnitTest/test_coarseSunSensor.py
index c4c93b8aec..53f97d9a8b 100644
--- a/src/simulation/sensors/coarseSunSensor/_UnitTest/test_coarseSunSensor.py
+++ b/src/simulation/sensors/coarseSunSensor/_UnitTest/test_coarseSunSensor.py
@@ -262,8 +262,8 @@ def setupCSS(CSS):
         constellationP1data = dataLogP1.CosValue
         constellationP2data = dataLogP2.CosValue
 
-        plt.figure(1, figsize=(7, 5), dpi=80, facecolor='w', edgecolor='k')
-        plt.clf()
+        if name == 'plain':
+            plt.figure(figsize=(7, 5), dpi=80, facecolor='w', edgecolor='k')
         plt.subplot(2, 1, 1)
         for i in range(4):
             sensorlabel = "cssP1" + str(i + 1)
@@ -293,7 +293,8 @@ def setupCSS(CSS):
     else:
         justTheNoise = cssOutput - truthVector  # subtract curve from noisy curve
         outputStd = np.std(justTheNoise)
-        plt.figure(3, figsize=(7, 5), dpi=80, facecolor='w', edgecolor='k')
+        if name == 'plain':
+            plt.figure(3, figsize=(7, 5), dpi=80, facecolor='w', edgecolor='k')
         plt.plot(dataLogSingle.times() * macros.NANO2MIN, cssOutput, label=name, zorder=zLevel, linewidth=lineWide)
         plt.legend()
         plt.xlabel('Time [min]')
diff --git a/src/simulation/sensors/imuSensor/_UnitTest/test_imu_sensor.py b/src/simulation/sensors/imuSensor/_UnitTest/test_imu_sensor.py
index e071d02b78..7a0d9a22d0 100755
--- a/src/simulation/sensors/imuSensor/_UnitTest/test_imu_sensor.py
+++ b/src/simulation/sensors/imuSensor/_UnitTest/test_imu_sensor.py
@@ -340,6 +340,7 @@ def unitSimIMU(show_plots,   testCase,       stopTime,       procRate, gyroLSBIn
     rDotDotOut = dataLog.AccelPlatform
     DVout = dataLog.DVFramePlatform
 
+    plt.close('all')
     # truth/output comparison plots and AutoTex output
     time = dataLog.times()/1e9
     plt.figure(1,figsize=(7, 5), dpi=80, facecolor='w', edgecolor='k')
@@ -522,7 +523,7 @@ def unitSimIMU(show_plots,   testCase,       stopTime,       procRate, gyroLSBIn
                                          'height=0.7\\textwidth, keepaspectratio', path)
         if show_plots:
             plt.show()
-            plt.close('all')
+        plt.close('all')
 
     #
     # Outputs to AutoTex
diff --git a/src/simulation/sensors/simpleVoltEstimator/_UnitTest/test_simpleVoltEstimator.py b/src/simulation/sensors/simpleVoltEstimator/_UnitTest/test_simpleVoltEstimator.py
index b136dcb353..9de33e240f 100755
--- a/src/simulation/sensors/simpleVoltEstimator/_UnitTest/test_simpleVoltEstimator.py
+++ b/src/simulation/sensors/simpleVoltEstimator/_UnitTest/test_simpleVoltEstimator.py
@@ -124,8 +124,7 @@ def unitSimpleVoltEstimator(show_plots):
             testFailCount += 1
             testMessages.append("FAILED: Too few error counts - " + str(count))
 
-    plt.figure(1)
-    plt.clf()
+    plt.close('all')
     plt.figure(1, figsize=(7, 5), dpi=80, facecolor='w', edgecolor='k')
     plt.plot(dataVoltLog.times() * 1.0E-9, volt[:])
 
@@ -135,7 +134,7 @@ def unitSimpleVoltEstimator(show_plots):
                                      r'height=0.4\textwidth, keepaspectratio', path)
     if show_plots:
         plt.show()
-        plt.close('all')
+    plt.close('all')
 
     # check if BSK_ERROR is returned if pMatrix is wrong size
     pMatrixBad = [[0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.],
diff --git a/src/tests/test_scenarioIntegrators.py b/src/tests/test_scenarioIntegrators.py
index 016df0e6d3..a1956ab89e 100644
--- a/src/tests/test_scenarioIntegrators.py
+++ b/src/tests/test_scenarioIntegrators.py
@@ -55,7 +55,7 @@ def test_scenarioIntegrators(show_plots):
     testFailCount = 0  # zero unit test result counter
     testMessages = []  # create empty array to store test log messages
 
-    for integratorCase in ["rk4", "rkf45", "rkf78", "euler", "rk2"]:
+    for integratorCase in ["rk4", "rkf45", "rkf78", "euler", "rk2", "bogackiShampine"]:
 
         # each test method requires a single assert method to be called
         posData, figureList = scenarioIntegrators.run(show_plots, integratorCase)
@@ -105,6 +105,14 @@ def test_scenarioIntegrators(show_plots):
                 ,[ 4614818.45321768,   -3602456.431072683,  -3837076.4216056713]
                 ,[ 5879286.370365726,   3561242.507948514,  -1319786.6261035257]
             ]
+        if integratorCase == "bogackiShampine":
+            truePos = [
+                [-2816801.601023492,  5248174.846916147,   3677157.2646772973]
+                , [-6380593.356301512,  -1467889.8216560283,  2481802.3429243616]
+                , [-2234294.4280201513, -6411056.176926798,  -1712849.5540222875]
+                , [ 4611359.466030476,  -3607692.026362182,  -3837674.903973699 ]
+                , [ 5882019.471581395,   3556295.378996462,  -1323290.9895207654]
+            ]
 
     # compare the results to the truth values
     accuracy = 1.0  # meters