diff --git a/pyoptsparse/pySLSQP/pySLSQP.py b/pyoptsparse/pySLSQP/pySLSQP.py
index c8da3852..0a13b0d7 100644
--- a/pyoptsparse/pySLSQP/pySLSQP.py
+++ b/pyoptsparse/pySLSQP/pySLSQP.py
@@ -220,6 +220,10 @@ def slgrad(m, me, la, n, f, g, df, dg, x):
             # fmt: on
             optTime = time.time() - t0
 
+            # Clip final result to user bounds (this occurs during the optimization as well
+            # so this just makes the output consistent with what the optimizer sees)
+            xs = np.clip(xs, blx, bux)
+
             # some entries of W include the lagrange multipliers
             # for each constraint, there are two entries (lower, upper).
             # if only one is active, look for the nonzero. If both are active, take the first one
diff --git a/tests/test_slsqp.py b/tests/test_slsqp.py
new file mode 100644
index 00000000..e8911466
--- /dev/null
+++ b/tests/test_slsqp.py
@@ -0,0 +1,47 @@
+"""Test class for SLSQP specific tests"""
+
+# First party modules
+from pyoptsparse import Optimization, OPT
+
+# Local modules
+from testing_utils import OptTest
+
+
+class TestSLSQP(OptTest):
+
+    def test_slsqp_strong_bound_enforcement(self):
+        """
+        Test that SLSQP will never evaluate the function or gradient outside
+        the design variable bounds. Without strong bound enforcement, the
+        optimizer will step outside the bounds and a ValueError will be raised.
+        With strong bound enforement, this code will run without raising any
+        errors
+        """
+
+        def objfunc(xdict):
+            x = xdict["xvars"]
+            funcs = {}
+            if x[0] < 0:
+                raise ValueError("Function cannot be evaluated below 0.")
+            funcs["obj"] = (x[0] - 1.0) ** 2
+            fail = False
+            return funcs, fail
+
+        def sens(xdict, funcs):
+            x = xdict["xvars"]
+            if x[0] < 0:
+                raise ValueError("Function cannot be evaluated below 0.")
+            funcsSens = {
+                "obj": {"xvars": [2 * (x[0] + 1.0)]},
+            }
+            fail = False
+            return funcsSens, fail
+
+        optProb = Optimization("Problem with Error Region", objfunc)
+        optProb.addVarGroup("xvars", 1, lower=[0], value=[2])
+        optProb.addObj("obj")
+        opt = OPT("SLSQP")
+        sol = opt(optProb, sens=sens)
+        self.assertEqual(sol.optInform["value"], 0)
+        self.assertGreaterEqual(sol.xStar["xvars"][0], 0)
+        self.assertAlmostEqual(sol.xStar["xvars"][0], 0, places=9)