diff --git a/src/iminuit/minuit.py b/src/iminuit/minuit.py
index 63d6deca..399d5525 100644
--- a/src/iminuit/minuit.py
+++ b/src/iminuit/minuit.py
@@ -987,6 +987,7 @@ def scipy(
         hess: Any = None,
         hessp: Any = None,
         constraints: Iterable = None,
+        options: Optional[Dict[str, Any]] = None,
     ) -> "Minuit":
         """
         Minimize with SciPy algorithms.
@@ -1014,6 +1015,10 @@ def scipy(
             as the original fcn, see hess parameter for details. No parameters may be
             omitted in the signature, even if those parameters are not used in the
             constraint.
+        options : dict, optional
+            A dictionary of solver options to pass to the SciPy minimizer through the
+            `options` parameter of :func:`scipy.optimize.minimize`. See each solver
+            method for the options it accepts.
 
         Notes
         -----
@@ -1027,6 +1032,13 @@ def scipy(
         criterion is evaluated only after the original algorithm already stopped. This
         means that usually SciPy minimizers will use more iterations than Migrad and
         the tolerance :attr:`tol` has no effect on SciPy minimizers.
+
+        You can specify convergence tolerance and other options for the SciPy minimizers
+        through the `options` parameter. Note that providing the SciPy options
+        `"maxiter"`, `"maxfev"`, and/or `"maxfun"` (depending on the minimizer) takes
+        precedence over providing a value for `ncall`. If you want to explicitly control
+        the number of iterations or function evaluations for a particular SciPy minimizer,
+        you should provide values for all of its relevant options.
         """
         try:
             from scipy.optimize import (
@@ -1224,18 +1236,30 @@ def __call__(self, par, v):
             else:
                 method = "BFGS"
 
+        options = options or {}
+
+        # attempt to set default number of function evaluations if not provided
         # various workarounds for API inconsistencies in scipy.optimize.minimize
-        options = {"maxiter": ncall}
+        added_maxiter = False
+        if "maxiter" not in options:
+            options["maxiter"] = ncall
+            added_maxiter = True
         if method in (
             "Nelder-Mead",
             "Powell",
         ):
-            options["maxfev"] = ncall
-            del options["maxiter"]
+            if "maxfev" not in options:
+                options["maxfev"] = ncall
+
+            if added_maxiter:
+                del options["maxiter"]
 
         if method in ("L-BFGS-B", "TNC"):
-            options["maxfun"] = ncall
-            del options["maxiter"]
+            if "maxfun" not in options:
+                options["maxfun"] = ncall
+
+            if added_maxiter:
+                del options["maxiter"]
 
         if method in ("COBYLA", "SLSQP", "trust-constr") and constraints is None:
             constraints = ()
diff --git a/tests/test_scipy.py b/tests/test_scipy.py
index c4c9957a..844ee949 100644
--- a/tests/test_scipy.py
+++ b/tests/test_scipy.py
@@ -1,6 +1,6 @@
 import pytest
 from numpy.testing import assert_allclose
-from iminuit import Minuit
+from iminuit import Minuit, cost
 from iminuit.testing import rosenbrock, rosenbrock_grad
 import numpy as np
 
@@ -252,3 +252,31 @@ def test_on_modified_state():
     m.scipy()  # used to fail
     assert m.valid
     assert_allclose(m.values, [0, 2], atol=1e-3)
+
+
+def test_options():
+    # simple example of uniform pdf with bounds on b to show tolerance
+    # can be improved with options
+    def density(x, b):
+        return b, np.full_like(x, b)
+
+    # with empty data, b=0
+    c = cost.ExtendedUnbinnedNLL([], density)
+
+    # Minimize with scipy's Powell and store the value of b
+    m = Minuit(c, b=0)
+    m.limits["b"] = (0, None)
+
+    m.scipy(method="Powell")
+    b_without_options = m.values["b"]
+
+    # try using scipy options to show it is better
+    c = cost.ExtendedUnbinnedNLL([], density)
+
+    m = Minuit(c, b=0)
+    m.limits["b"] = (0, None)
+
+    m.scipy(method="Powell", options={"xtol": 1e-10, "ftol": 1e-10})
+    b_with_options = m.values["b"]
+
+    assert b_without_options > b_with_options