diff --git a/benchmarks/femc_electron/analysis/femc_electron_plots.py b/benchmarks/femc_electron/analysis/femc_electron_plots.py
index cc34633..3a7128b 100644
--- a/benchmarks/femc_electron/analysis/femc_electron_plots.py
+++ b/benchmarks/femc_electron/analysis/femc_electron_plots.py
@@ -132,7 +132,7 @@ def gauss(x, A,mu, sigma):
     sigma=np.sqrt(y[slc])+0.5*(y[slc]==0)
     p0=(100, p, 3)
 
-    coeff, var_matrix = curve_fit(fnc, list(bcs[slc]), list(y[slc]), p0=p0,sigma=list(sigma))
+    coeff, var_matrix = curve_fit(fnc, list(bcs[slc]), list(y[slc]), p0=p0,sigma=list(sigma), maxfev=10000)
     #res=np.abs(coeff[2]/coeff[1])
     if p==50:
         xx=np.linspace(15*p/20,22*p/20, 100)
@@ -151,7 +151,7 @@ def gauss(x, A,mu, sigma):
 plt.errorbar(pvals, 100*np.array(res), 100*np.array(dres), ls='', marker='o')
 fnc = lambda E, a, b: np.hypot(a,b/np.sqrt(E))
 p0=(.05, .12)
-coeff, var_matrix = curve_fit(fnc, pvals, res, p0=p0,sigma=dres)
+coeff, var_matrix = curve_fit(fnc, pvals, res, p0=p0,sigma=dres, maxfev=10000)
 xx=np.linspace(7, 85, 100)
 plt.plot(xx, 100*fnc(xx,*coeff), label=f'fit:{100*coeff[0]:.1f}%$\\oplus\\frac{{{100*coeff[1]:.0f}\\%}}{{\\sqrt{{E}}}}$')
 plt.legend()
@@ -190,7 +190,7 @@ def gauss(x, A,mu, sigma):
         p0=(100, p, 3)
 
         try:
-            coeff, var_matrix = curve_fit(fnc, list(bcs[slc]), list(y[slc]), p0=p0,sigma=list(sigma))
+            coeff, var_matrix = curve_fit(fnc, list(bcs[slc]), list(y[slc]), p0=p0,sigma=list(sigma), maxfev=10000)
             if abs(coeff[1])>100 or np.sqrt(var_matrix[1][1])>100:
                 continue
             pvals.append(p)
diff --git a/benchmarks/femc_photon/analysis/femc_photon_plots.py b/benchmarks/femc_photon/analysis/femc_photon_plots.py
index 8170f4a..f08ca2d 100644
--- a/benchmarks/femc_photon/analysis/femc_photon_plots.py
+++ b/benchmarks/femc_photon/analysis/femc_photon_plots.py
@@ -131,7 +131,7 @@ def gauss(x, A,mu, sigma):
     sigma=np.sqrt(y[slc])+0.5*(y[slc]==0)
     p0=(100, p, 3)
 
-    coeff, var_matrix = curve_fit(fnc, list(bcs[slc]), list(y[slc]), p0=p0,sigma=list(sigma))
+    coeff, var_matrix = curve_fit(fnc, list(bcs[slc]), list(y[slc]), p0=p0,sigma=list(sigma), maxfev=10000)
     #res=np.abs(coeff[2]/coeff[1])
     if p==50:
         xx=np.linspace(15*p/20,22*p/20, 100)
@@ -150,7 +150,7 @@ def gauss(x, A,mu, sigma):
 plt.errorbar(pvals, 100*np.array(res), 100*np.array(dres), ls='', marker='o')
 fnc = lambda E, a, b: np.hypot(a,b/np.sqrt(E))
 p0=(.05, .12)
-coeff, var_matrix = curve_fit(fnc, pvals, res, p0=p0,sigma=dres)
+coeff, var_matrix = curve_fit(fnc, pvals, res, p0=p0,sigma=dres, maxfev=10000)
 xx=np.linspace(7, 85, 100)
 plt.plot(xx, 100*fnc(xx,*coeff), label=f'fit:{100*coeff[0]:.1f}%$\\oplus\\frac{{{100*coeff[1]:.0f}\\%}}{{\\sqrt{{E}}}}$')
 plt.legend()
@@ -188,7 +188,7 @@ def gauss(x, A,mu, sigma):
         sigma=np.sqrt(y[slc])+0.5*(y[slc]==0)
         p0=(100, p, 3)
         try:
-            coeff, var_matrix = curve_fit(fnc, list(bcs[slc]), list(y[slc]), p0=p0,sigma=list(sigma))
+            coeff, var_matrix = curve_fit(fnc, list(bcs[slc]), list(y[slc]), p0=p0,sigma=list(sigma), maxfev=10000)
             if abs(coeff[1])>100 or np.sqrt(var_matrix[1][1])>100:
                 continue
             pvals.append(p)
diff --git a/benchmarks/femc_pi0/analysis/femc_pi0_plots.py b/benchmarks/femc_pi0/analysis/femc_pi0_plots.py
index cdf9d23..a732c31 100644
--- a/benchmarks/femc_pi0/analysis/femc_pi0_plots.py
+++ b/benchmarks/femc_pi0/analysis/femc_pi0_plots.py
@@ -167,7 +167,7 @@ def gauss(x, A,mu, sigma):
     sigma=np.sqrt(y[slc])+0.5*(y[slc]==0)
     p0=(100, p, 3)
 
-    coeff, var_matrix = curve_fit(fnc, list(bcs[slc]), list(y[slc]), p0=p0,sigma=list(sigma))
+    coeff, var_matrix = curve_fit(fnc, list(bcs[slc]), list(y[slc]), p0=p0,sigma=list(sigma), maxfev=10000)
     #res=np.abs(coeff[2]/coeff[1])
     if p==50:
         xx=np.linspace(15*p/20,22*p/20, 100)
@@ -186,7 +186,7 @@ def gauss(x, A,mu, sigma):
 plt.errorbar(pvals, 100*np.array(res), 100*np.array(dres), ls='', marker='o')
 fnc = lambda E, a, b: np.hypot(a,b/np.sqrt(E))
 p0=(.05, .12)
-coeff, var_matrix = curve_fit(fnc, pvals, res, p0=p0,sigma=dres)
+coeff, var_matrix = curve_fit(fnc, pvals, res, p0=p0,sigma=dres, maxfev=10000)
 xx=np.linspace(15, 85, 100)
 plt.plot(xx, 100*fnc(xx,*coeff), label=f'fit:{100*coeff[0]:.1f}%$\\oplus\\frac{{{100*coeff[1]:.0f}\\%}}{{\\sqrt{{E}}}}$')
 plt.legend()
@@ -225,7 +225,7 @@ def gauss(x, A,mu, sigma):
         p0=(100, p, 3)
 
         try:
-            coeff, var_matrix = curve_fit(fnc, list(bcs[slc]), list(y[slc]), p0=p0,sigma=list(sigma))
+            coeff, var_matrix = curve_fit(fnc, list(bcs[slc]), list(y[slc]), p0=p0,sigma=list(sigma), maxfev=10000)
             if abs(coeff[1])>100 or np.sqrt(var_matrix[1][1])>100:
                 continue
             pvals.append(p)