diff --git a/cyipopt/scipy_interface.py b/cyipopt/scipy_interface.py
index a601ca96..2bbe4fb8 100644
--- a/cyipopt/scipy_interface.py
+++ b/cyipopt/scipy_interface.py
@@ -26,6 +26,7 @@
     SCIPY_INSTALLED = True
     del scipy
     from scipy.optimize import approx_fprime
+    from scipy.optimize._numdiff import approx_derivative
     try:
         from scipy.optimize import OptimizeResult
     except ImportError:
@@ -99,7 +100,7 @@ def __init__(self, fun, args=(), kwargs=None, jac=None, hess=None,
             con_args = con.get('args', [])
             con_kwargs = con.get('kwargs', [])
             if con_jac is None:
-                con_jac = lambda x0, *args, **kwargs: approx_fprime(x0, con_fun, eps, *args, **kwargs)
+                con_jac = lambda x: approx_derivative(con_fun, x, method='3-point')
             self._constraint_funs.append(con_fun)
             self._constraint_jacs.append(con_jac)
             self._constraint_args.append(con_args)
diff --git a/cyipopt/tests/unit/test_scipy_optional.py b/cyipopt/tests/unit/test_scipy_optional.py
index ce5dd245..431fc85b 100644
--- a/cyipopt/tests/unit/test_scipy_optional.py
+++ b/cyipopt/tests/unit/test_scipy_optional.py
@@ -15,6 +15,7 @@
 import cyipopt
 
 
+
 @pytest.mark.skipif("scipy" in sys.modules,
                     reason="Test only valid if no Scipy available.")
 def test_minimize_ipopt_import_error_if_no_scipy():
@@ -59,14 +60,32 @@ def test_minimize_ipopt_nojac_if_scipy():
                     reason="Test only valid if Scipy available.")
 def test_minimize_ipopt_nojac_constraints_if_scipy():
     """ `minimize_ipopt` works without Jacobian and with constraints"""
-    from scipy.optimize import rosen
+    from scipy.optimize import rosen, rosen_der
     x0 = [1.3, 0.7, 0.8, 1.9, 1.2]
     constr = {"fun": lambda x: rosen(x) - 1.0, "type": "ineq"}
-    res = cyipopt.minimize_ipopt(rosen, x0, constraints=constr)
+    res = cyipopt.minimize_ipopt(rosen, x0, jac=rosen_der, constraints=constr)
     assert isinstance(res, dict)
     assert np.isclose(res.get("fun"), 1.0)
     assert res.get("status") == 0
     assert res.get("success") is True
-    expected_res = np.array([1.001867, 0.99434067, 1.05070075, 1.17906312,
-                             1.38103001])
+    expected_res = np.array([1.00407015, 0.99655763, 1.05556205, 1.18568342, 1.38386505])
+    np.testing.assert_allclose(res.get("x"), expected_res)
+
+
+@pytest.mark.skipif("scipy" not in sys.modules,
+                    reason="Test only valid if Scipy available.")
+def test_minimize_ipopt_nojac_constraints_vectorvalued_if_scipy():
+    """ `minimize_ipopt` works without Jacobian and with constraints 
+    without specifying the jacobian of the vector valued constraint 
+    function"""
+    from scipy.optimize import rosen, rosen_der
+    x0 = np.array([0.5, 0.75])
+    bounds = [np.array([0, 1]), np.array([-0.5, 2.0])]
+    expected_res = 0.25 * np.ones_like(x0)
+    eq_cons = {"fun" : lambda x: x - expected_res, "type": "eq"}
+    res = cyipopt.minimize_ipopt(rosen, x0, jac=rosen_der, bounds=bounds, constraints=[eq_cons])
+    assert isinstance(res, dict)
+    assert res.get("status") == 0
+    assert res.get("success") is True
     np.testing.assert_allclose(res.get("x"), expected_res)
+    
\ No newline at end of file
diff --git a/cyipopt/utils.py b/cyipopt/utils.py
index 048ae2d1..822c8577 100644
--- a/cyipopt/utils.py
+++ b/cyipopt/utils.py
@@ -7,8 +7,6 @@
 import warnings
 from functools import wraps
 
-from ipopt_wrapper import Problem
-
 
 def deprecated_warning(new_name):
     """Decorator that issues a FutureWarning for deprecated functionality.
@@ -91,3 +89,4 @@ def generate_deprecation_warning_msg(what,
            f"deprecated in CyIpopt. Please replace all uses and use "
            f"'{new_name}' going forward.")
     return msg
+
diff --git a/examples/constraints.py b/examples/constraints.py
new file mode 100644
index 00000000..32e159ee
--- /dev/null
+++ b/examples/constraints.py
@@ -0,0 +1,14 @@
+import numpy as np
+from scipy.optimize import rosen, rosen_der
+from cyipopt import minimize_ipopt
+
+x0 = np.array([0.5, 0])
+
+bounds = [np.array([0, 1]), np.array([-0.5, 2.0])]
+
+eq_cons = {"type": "eq",
+           "fun": lambda x: np.array([2*x[0] + x[1] - 1, x[0]**2 - 0.1])}
+
+res = minimize_ipopt(rosen, x0, jac=rosen_der, bounds=bounds, constraints=[eq_cons])
+
+print(res)