diff --git a/deerlab/classes.py b/deerlab/classes.py
index 39063cbf..49a6ec46 100644
--- a/deerlab/classes.py
+++ b/deerlab/classes.py
@@ -522,7 +522,7 @@ def propagate(self,model,lb=None,ub=None,samples=None):
         ubm : ndarray
             Upper bounds of the values returned by ``model``, by default assumed unconstrained.
         samples : int, optional
-            Number of samples to use when propagating uncertainty. If not provided, default value is 1000.
+            Number of samples to use when propagating a bootstraped uncertainty. If not provided, default value is 1000.
         
         Returns
         -------
@@ -583,7 +583,7 @@ def propagate(self,model,lb=None,ub=None,samples=None):
                 # Get the parameter uncertainty distribution
                 values,pdf = self.pardist(n)
                 # Random sampling form the uncertainty distribution
-                sampled_parameters[n] =  [np.random.choice(values, p=pdf/sum(pdf)) for _ in range(Nsamples)]
+                sampled_parameters[n] = np.random.choice(values, p=pdf/sum(pdf),size=Nsamples)
             # Convert to matrix
             sampled_parameters = np.atleast_2d(sampled_parameters)
 
diff --git a/deerlab/dd_models.py b/deerlab/dd_models.py
index 1c4d6d3f..e590f7c2 100644
--- a/deerlab/dd_models.py
+++ b/deerlab/dd_models.py
@@ -99,7 +99,7 @@ def docstr_example(fcnstr):
 # =================================================================
 def _normalize(r,P):
     if not all(P==0):
-        P = P/np.trapz(P,r)
+        P = P/np.trapezoid(P,r)
     return P
 # =================================================================
 
@@ -110,7 +110,7 @@ def _multigaussfun(r,r0,sig):
     P = np.sqrt(1/(2*np.pi))*1/sig*np.exp(-0.5*((r.T-r0)/sig)**2)
     if not np.all(P==0):
         # Normalization
-        P = np.squeeze(P)/np.sum([np.trapezoid(c,r) for c in P.T])
+        P = np.squeeze(P)/np.array([np.trapezoid(c,r) for c in P.T]).sum()
     else: 
         P = np.squeeze(P)
     return P
diff --git a/deerlab/fit.py b/deerlab/fit.py
index 7c757649..71d26231 100644
--- a/deerlab/fit.py
+++ b/deerlab/fit.py
@@ -361,7 +361,9 @@ def fit(model_, y, *constants, par0=None, penalties=None, bootstrap=0, noiselvl=
     paramlist : list
         List of the fitted parameter names ordered according to the model parameter indices.
     model : ndarray
-        Fitted model response.     
+        Fitted model response.   
+    modelUncert : :ref:`UQResult`
+        Uncertainty quantification of the fitted model response.
     regparam : scalar
         Regularization parameter value used for the regularization of the linear parameters.
     penweights : scalar or list thereof 
@@ -470,7 +472,7 @@ def fit(model_, y, *constants, par0=None, penalties=None, bootstrap=0, noiselvl=
     fitfcn = lambda y,penweights: snlls(y, Amodel_fcn, par0, lb=lb, ub=ub, lbl=lbl, ubl=ubl, mask=mask, weights=weights, 
                                                 subsets=ysubsets, lin_frozen=linfrozen, nonlin_frozen=nonlinfrozen,
                                                 regparam=regparam, reg=reg, regparamrange=regparamrange, noiselvl=noiselvl,
-                                                extrapenalty=extrapenalties(penweights), **kwargs)        
+                                                extrapenalty=extrapenalties(penweights), modeluq=True, **kwargs)        
 
     # Prepare outer optimization of the penalty weights, if necessary
     fitfcn = _outerOptimization(fitfcn,penalties,sigmas)
@@ -491,14 +493,19 @@ def bootstrap_fcn(ysim):
         else:
             bootstrap_verbose = False
             
-        param_uq = bootstrap_analysis(bootstrap_fcn,ysplit,fitresults.model,samples=bootstrap,noiselvl=noiselvl,cores=bootcores, verbose=bootstrap_verbose)
+        param_uq = bootstrap_analysis(bootstrap_fcn,ysplit,fitresults.model,samples=bootstrap-1,noiselvl=noiselvl,cores=bootcores, verbose=bootstrap_verbose)
         # Include information on the boundaries for better uncertainty estimates
         paramlb = model._vecsort(model._getvector('lb'))[np.concatenate(param_idx)] 
         paramub = model._vecsort(model._getvector('ub'))[np.concatenate(param_idx)] 
         fitresults.paramUncert = UQResult('bootstrap',data=param_uq[0].samples,lb=paramlb,ub=paramub)
         fitresults.param = fitresults.paramUncert.median
+
         # Get the uncertainty estimates for the model response
+        modellb = np.min(param_uq[1].samples,axis=0)
+        modelub = np.max(param_uq[1].samples,axis=0)
+            
         fitresults.model = [param_uq[n].median for n in range(1,len(param_uq))]
+        fitresults.modelUncert = UQResult('bootstrap',data=param_uq[1].samples,lb=modellb,ub=modelub)
         if len(fitresults.model)==1: 
             fitresults.model = fitresults.model[0]
     # Get some basic information on the parameter vector
@@ -509,7 +516,7 @@ def bootstrap_fcn(ysim):
     # Dictionary of parameter names and fit uncertainties
     FitResult_paramuq = {f'{key}Uncert': model._getparamuq(fitresults.paramUncert,idx) for key,idx in zip(keys,param_idx)}
     # Dictionary of other fit quantities of interest
-    FitResult_dict = {key: getattr(fitresults,key) for key in ['y','mask','param','paramUncert','model','cost','plot','residuals','stats','regparam','regparam_stats','__plot_inputs']}
+    FitResult_dict = {key: getattr(fitresults,key) for key in ['y','mask','param','paramUncert','model','modelUncert','cost','plot','residuals','stats','regparam','regparam_stats','__plot_inputs']}
     _paramlist = model._parameter_list('vector')
 
     param_idx = [[] for _ in _paramlist]
@@ -536,8 +543,8 @@ def _scale(x):
                     FitResult_param_[f'{key}_scale'] = _scale(FitResult_param_[key]) # Normalization factor
                     FitResult_param_[key] = param.normalization(FitResult_param_[key]) # Normalized value
 
-                    FitResult_paramuq_[f'{key}_scaleUncert'] = FitResult_paramuq_[f'{key}Uncert'].propagate(_scale)
-                    FitResult_paramuq_[f'{key}Uncert'] = FitResult_paramuq_[f'{key}Uncert'].propagate(lambda x: x/FitResult_param_[f'{key}_scale'], lb=param.lb, ub=param.ub) # Normalization of the uncertainty
+                    FitResult_paramuq_[f'{key}_scaleUncert'] = FitResult_paramuq_[f'{key}Uncert'].propagate(_scale,samples=bootstrap)
+                    FitResult_paramuq_[f'{key}Uncert'] = FitResult_paramuq_[f'{key}Uncert'].propagate(lambda x: x/FitResult_param_[f'{key}_scale'], lb=param.lb, ub=param.ub,samples=bootstrap) # Normalization of the uncertainty
     if len(noiselvl)==1: 
         noiselvl = noiselvl[0]
     
diff --git a/deerlab/fitresult.py b/deerlab/fitresult.py
index 4841f9e6..27d33aa9 100644
--- a/deerlab/fitresult.py
+++ b/deerlab/fitresult.py
@@ -97,20 +97,21 @@ def _extarct_params_from_model(self, model):
         if not hasattr(self,'param'):
             raise ValueError('The fit object does not contain any fitted parameters.')
 
-        # # Enforce model normalization
-        # normfactor_keys = []
-        # for key in modelparam:
-        #     param = getattr(model,key)
-        #     if np.all(param.linear):
-        #         if param.normalization is not None:
-        #             normfactor_key = f'{key}_scale'
-        #             normfactor_keys.append(normfactor_key)
-        #             try:
-        #                 model.addnonlinear(normfactor_key,lb=-np.inf,ub=np.inf,par0=1,description=f'Normalization factor of {key}')
-        #                 getattr(model,normfactor_key).freeze(1)
-        #             except KeyError:
-        #                 pass
-                    
+        # Enforce model normalization
+        normfactor_keys = []
+        for key in modelparam:
+            param = getattr(model,key)
+            if np.all(param.linear):
+                if param.normalization is not None:
+                    normfactor_key = f'{key}_scale'
+                    normfactor_keys.append(normfactor_key)
+                    try:
+                        model.addnonlinear(normfactor_key,lb=-np.inf,ub=np.inf,par0=1,description=f'Normalization factor of {key}')
+                        getattr(model,normfactor_key).freeze(1)
+                    except KeyError:
+                        pass
+        modelparam += normfactor_keys
+        
 
         # # Get some basic information on the parameter vector
         # modelparam = model._parameter_list(order='vector')
@@ -187,6 +188,7 @@ def evaluate(self, model, *constants):
             Model response at the fitted parameter values. 
         """
         try:
+            model = model.copy()
             modelparam = model._parameter_list('vector')
             modelparam, fitparams, fitparam_idx = self._extarct_params_from_model(model)
         except AttributeError:
@@ -199,7 +201,7 @@ def evaluate(self, model, *constants):
         response = model(*constants,**parameters)
         return response
     
-    def propagate(self, model, *constants, lb=None, ub=None):
+    def propagate(self, model, *constants, lb=None, ub=None,samples=None):
         """
         Propagate the uncertainty in the fit results to a model's response.
 
@@ -223,7 +225,10 @@ def propagate(self, model, *constants, lb=None, ub=None):
         lb : array_like, optional 
             Lower bounds of the model response.
         ub : array_like, optional 
-            Upper bounds of the model response.   
+            Upper bounds of the model response. 
+        samples : int, optional
+            Number of samples to use when propagating a bootstraped uncertainty. If not provided, default value is 1000.
+
 
         Returns
         -------
@@ -231,8 +236,10 @@ def propagate(self, model, *constants, lb=None, ub=None):
         responseUncert : :ref:`UQResult`
             Uncertainty quantification of the model's response.
         """
+        
 
         try:
+            model = model.copy()
             modelparam = model._parameter_list('vector')
             modelparam, fitparams, fitparam_idx = self._extarct_params_from_model(model)
 
@@ -241,7 +248,7 @@ def propagate(self, model, *constants, lb=None, ub=None):
 
 
         # Propagate the uncertainty from that subset to the model
-        modeluq = self.paramUncert.propagate(lambda param: model(*constants,*[param[s] for s in fitparam_idx]),lb,ub)
+        modeluq = self.paramUncert.propagate(lambda param: model(*constants,*[param[s] for s in fitparam_idx]),lb,ub,samples)
         return modeluq
     
 
diff --git a/deerlab/model.py b/deerlab/model.py
index 7ba709dd..5e4866a3 100644
--- a/deerlab/model.py
+++ b/deerlab/model.py
@@ -218,6 +218,13 @@ def unfreeze(self):
             self.frozen = False
             self.value = None
     #---------------------------------------------------------------------------------------
+    def copy(self):
+        """
+        Return a deep copy of the parameter
+        """
+
+        return deepcopy(self)
+
 #===================================================================================
 
 
@@ -867,7 +874,7 @@ def __call__(self,*args,**kargs):
 
         # Check that all parameters have been passed
         if len(θ)!=self.Nparam:
-            raise SyntaxError(f'The model requires {self.Nparam} parameters, but {len(args_list)} were specified.')   
+            raise SyntaxError(f'The model requires {self.Nparam} parameters, but {len(θ)} were specified.')   
 
         # Determine which parameters are linear and which nonlinear
         θlin, θnonlin = self._split_linear(θ)
@@ -973,6 +980,13 @@ def __str__(self):
         """
         return self._parameter_table()      
     #---------------------------------------------------------------------------------------
+    def copy(self):
+        """
+        Return a deep copy of the model
+        """
+
+        return deepcopy(self)
+
 #===================================================================================
 
 #==============================================================================
diff --git a/examples/advanced/ex_dipolarpathways_selection.py b/examples/advanced/ex_dipolarpathways_selection.py
index 57fb0e18..bf255174 100644
--- a/examples/advanced/ex_dipolarpathways_selection.py
+++ b/examples/advanced/ex_dipolarpathways_selection.py
@@ -11,7 +11,6 @@
 import numpy as np
 import matplotlib.pyplot as plt
 from matplotlib.gridspec import GridSpec
-from mpl_toolkits.axes_grid1.inset_locator import inset_axes, InsetPosition
 import deerlab as dl
 
 
@@ -79,9 +78,7 @@
     Pfit = fits[n].P
     Pci = fits[n].PUncert.ci(95)
     # Setup the inset plot
-    axins = inset_axes(ax1,width="30%", height="30%", loc='upper left')
-    ip = InsetPosition(ax1,[0.35, 0.17+0.24*n, 0.6, 0.1])
-    axins.set_axes_locator(ip)
+    axins = ax1.inset_axes([0.35, 0.17+0.24*n, 0.6, 0.1])
     axins.yaxis.set_ticklabels([])
     axins.yaxis.set_visible(False)
     # Plot the distance distributions and their confidence bands
diff --git a/examples/advanced/ex_profileanalysis.py b/examples/advanced/ex_profileanalysis.py
index 82d5fc42..6a195bd5 100644
--- a/examples/advanced/ex_profileanalysis.py
+++ b/examples/advanced/ex_profileanalysis.py
@@ -48,7 +48,7 @@
 #%%
 # Extract fitted dipolar signal
 Vfit = results.model
-Vci = results.propagate(Vmodel).ci(95)
+Vci = results.modelUncert.ci(95)
 
 # Extract fitted distance distribution
 Pfit = results.P
diff --git a/examples/basic/ex_bootstrapping.py b/examples/basic/ex_bootstrapping.py
index de9a558f..196bd772 100644
--- a/examples/basic/ex_bootstrapping.py
+++ b/examples/basic/ex_bootstrapping.py
@@ -59,7 +59,7 @@
 
 # Extract fitted dipolar signal
 Vfit = results.model
-Vci = results.propagate(Vmodel).ci(95)
+Vci = results.modelUncert.ci(95)
 
 # Extract fitted distance distribution
 Pfit = results.P
@@ -78,6 +78,7 @@
 plt.plot(t,Vexp,'.',color='grey',label='Data')
 # Plot the fitted signal 
 plt.plot(t,Vfit,linewidth=3,label='Bootstrap median',color=violet)
+plt.fill_between(t,Vci[:,0],Vci[:,1],linewidth=0.1,label='Bootstrap median',color=violet,alpha=0.3)
 plt.plot(t,Bfit,'--',linewidth=3,color=violet,label='Unmodulated contribution')
 plt.legend(frameon=False,loc='best')
 plt.xlabel('Time $t$ (μs)')
diff --git a/examples/basic/ex_fitting_5pdeer.py b/examples/basic/ex_fitting_5pdeer.py
index 83d7519d..432d3129 100644
--- a/examples/basic/ex_fitting_5pdeer.py
+++ b/examples/basic/ex_fitting_5pdeer.py
@@ -58,7 +58,7 @@
 
 # Extract fitted dipolar signal
 Vfit = results.model
-Vci = results.propagate(Vmodel).ci(95)
+Vci = results.modelUncert.ci(95)
 
 # Extract fitted distance distribution
 Pfit = results.P
diff --git a/examples/intermediate/ex_fitting_5pdeer_pathways.py b/examples/intermediate/ex_fitting_5pdeer_pathways.py
index 370a7d93..516e50dc 100644
--- a/examples/intermediate/ex_fitting_5pdeer_pathways.py
+++ b/examples/intermediate/ex_fitting_5pdeer_pathways.py
@@ -47,7 +47,7 @@
 
 # Extract fitted dipolar signal
 Vfit = results.model
-Vci = results.propagate(Vmodel).ci(95)
+Vci = results.modelUncert.ci(95)
 
 # Extract fitted distance distribution
 Pfit = results.P
diff --git a/setup.py b/setup.py
index a0747967..5d237fc5 100644
--- a/setup.py
+++ b/setup.py
@@ -25,7 +25,7 @@
                         'joblib>=1.0.0',
                         'dill>=0.3.0',
                         'tqdm>=4.51.0',
-                        'matplotlib>=3.3.4',
+                        'matplotlib>=3.6.0',
                         'memoization>=0.3.1',
                         'pytest>=6.2.2',
                         'setuptools>=53.0.0',
diff --git a/test/test_model_class.py b/test/test_model_class.py
index b8f6b864..1a4c8a98 100644
--- a/test/test_model_class.py
+++ b/test/test_model_class.py
@@ -42,7 +42,14 @@ def mock_x():
 
 @pytest.fixture(scope='module')
 def mock_data(mock_x):
-    return bigauss(mock_x,mean1=3,mean2=4,std1=0.5,std2=0.2,amp1=0.5,amp2=0.6)
+    data = bigauss(mock_x,mean1=3,mean2=4,std1=0.5,std2=0.2,amp1=0.5,amp2=0.6)
+    return data
+
+@pytest.fixture(scope='module')
+def mock_data_noise(mock_x):
+    data = bigauss(mock_x,mean1=3,mean2=4,std1=0.5,std2=0.2,amp1=0.5,amp2=0.6)
+    data += whitegaussnoise(mock_x,0.01,seed=1)
+    return data
 
 model_types = ['parametric','semiparametric','semiparametric_vec',
                 'nonparametric','nonparametric_vec','nonparametric_vec_normalized']
@@ -511,7 +518,7 @@ def test_fit_model_confidence_intervals(mock_data,model_type,method):
     model = _generate_model(model_type, fixed_axis=True)
 
     if method=='bootstrap':
-        results = fit(model,mock_data + whitegaussnoise(x,0.01,seed=1), bootstrap=3)
+        results = fit(model,mock_data + whitegaussnoise(x,0.01,seed=1), bootstrap=100)
     else: 
         results = fit(model,mock_data + whitegaussnoise(x,0.01,seed=1))
     
@@ -546,6 +553,27 @@ def test_fit_evaluate_model(mock_data,mock_x,model_type):
         response *= results.scale
     
     assert np.allclose(response,mock_data)
+
+# ================================================================
+@pytest.mark.parametrize('method', ['bootstrap','moment'])
+@pytest.mark.parametrize('model_type', model_types)
+def test_fit_modelUncert(mock_data_noise,mock_x,model_type,method): 
+    "Check that the uncertainty of fit results can be calculated and is the uncertainty of the model is non zero for all but nonparametric models"
+    model = _generate_model(model_type, fixed_axis=False)
+    
+    if method=='bootstrap':
+        results = fit(model,mock_data_noise,mock_x, bootstrap=3)
+    else: 
+        results = fit(model,mock_data_noise,mock_x)
+
+    assert hasattr(results,'modelUncert')
+    ci_lower = results.modelUncert.ci(95)[:,0]
+    ci_upper = results.modelUncert.ci(95)[:,1]
+    assert np.less_equal(ci_lower,ci_upper).all()
+    if model_type != 'nonparametric' and model_type != 'nonparametric_vec' and model_type != 'nonparametric_vec_normalized':
+        assert np.all(np.round(ci_lower) <= np.round(results.model)) and np.less(ci_lower.sum(),results.model.sum())
+        assert np.all(np.round(ci_upper,5) >= np.round(results.model,5)) and np.greater(ci_upper.sum(),results.model.sum())
+
 # ================================================================
 
 # ================================================================