diff --git a/maoppy/psffit.py b/maoppy/psffit.py
index 2ac80834340a4508fc743e9416f28ee3dca7cdd7..3d0aaa26bcbe390fe158de366c12c1424b205160 100644
--- a/maoppy/psffit.py
+++ b/maoppy/psffit.py
@@ -49,8 +49,8 @@ def lsq_flux_bck(model, data, weights, background=True, positive_bck=False):
     return amp, bck
 
 
-def psffit(psf,model,x0,weights=None,dxdy=(0.,0.),flux_bck=(True,True),
-           positive_bck=False,fixed=None,npixfit=None,**kwargs):
+def psffit(psf, model, x0, weights=None, dxdy=(0.,0.), flux_bck=(True,True),
+           positive_bck=False, fixed=None, npixfit=None, max_nfev=None, **kwargs):
     """Fit a PSF with a parametric model solving the least-square problem
        epsilon(x) = SUM_pixel { weights * (amp * Model(x) + bck - psf)² }
     
@@ -77,6 +77,8 @@ def psffit(psf,model,x0,weights=None,dxdy=(0.,0.),flux_bck=(True,True),
         Fix some parameters to their initial value (default: None)
     npixfit : int (default=None)
         Increased pixel size for improved fitting accuracy
+    max_nfev : int or None
+        Maximum number of evaluation of the cost function
     **kwargs :
         All keywords used to call your `model`
     
@@ -172,10 +174,9 @@ def psffit(psf,model,x0,weights=None,dxdy=(0.,0.),flux_bck=(True,True),
         b_up = np.concatenate((b_up,[np.inf,np.inf]))
         return (b_low,b_up)
     
-    result = least_squares(cost, input2mini(x0,dxdy), bounds=get_bound(model))
+    result = least_squares(cost, input2mini(x0,dxdy), bounds=get_bound(model), max_nfev=max_nfev)
     
-    msg = result.message[0:min(len(result.message),15)]
-    print("\rPSFFIT finished %3u iter  %s..."%(cost.iter,msg))
+    print("\rPSFFIT finished in %u iter : %s"%(cost.iter,result.message))
     
     result.x, result.dxdy = mini2input(result.x) #split output between x and dxdy