diff --git a/gammapy/modeling/models/spectral.py b/gammapy/modeling/models/spectral.py
index 4e786b1cfd..7bc7512358 100644
--- a/gammapy/modeling/models/spectral.py
+++ b/gammapy/modeling/models/spectral.py
@@ -491,6 +491,8 @@ def plot_error(
         sed_type="dnde",
         energy_power=0,
         n_points=100,
+        n_samples=10000,
+        method="lin",
         **kwargs,
     ):
         """Plot spectral model error band.
@@ -536,6 +538,7 @@ def plot_error(
         ``n_points`` bins between the given bounds.
         """
         from gammapy.estimators.map.core import DEFAULT_UNIT
+        from gammapy.estimators import FluxMaps
 
         if self.is_norm_spectral_model:
             sed_type = "norm"
@@ -560,8 +563,56 @@ def plot_error(
             ax.yaxis.set_units(DEFAULT_UNIT[sed_type] * energy.unit**energy_power)
 
         flux, flux_err = self._get_plot_flux(sed_type=sed_type, energy=energy)
-        y_lo = scale_plot_flux(flux - flux_err, energy_power).quantity[:, 0, 0]
-        y_hi = scale_plot_flux(flux + flux_err, energy_power).quantity[:, 0, 0]
+
+        if method == "log":
+            flux_err_log = flux_err.data / flux.data
+            y_lo = flux / np.exp(flux_err_log)
+            y_hi = flux * np.exp(flux_err_log)
+            y_lo = scale_plot_flux(y_lo, energy_power).quantity[:, 0, 0]
+            y_hi = scale_plot_flux(y_hi, energy_power).quantity[:, 0, 0]
+        elif method == "lin":
+            y_lo = scale_plot_flux(flux - flux_err, energy_power).quantity[:, 0, 0]
+            y_hi = scale_plot_flux(flux + flux_err, energy_power).quantity[:, 0, 0]
+        elif method == "samples":
+            samples = np.random.multivariate_normal(
+                self.parameters.value, self.covariance, n_samples
+            )
+            # Reapply units to the sampled parameters
+            parameter_units = [par.unit for par in self.parameters]
+            samples = np.array([sample * parameter_units for sample in samples])
+
+            f = self(energy.center)
+            spl_ax = MapAxis(range(n_samples), node_type="center", name="sample")
+            f_samples = RegionNDMap.create(region=None, axes=[energy, spl_ax])
+            f_samples.quantity = (
+                np.array(
+                    [
+                        self.evaluate(energy.center, *samples[k, :])
+                        for k in range(n_samples)
+                    ]
+                )
+                * f.unit
+            )
+
+            f_samples_map = FluxMaps.from_maps(
+                dict(dnde=f_samples), sed_type="dnde", reference_model=self
+            )
+            f_samples = f_samples_map[sed_type]
+
+            y_lo = RegionNDMap.create(
+                region=None,
+                axes=[energy],
+                unit=f_samples.unit,
+                data=np.percentile(f_samples, 16, axis=0),
+            )
+            y_hi = RegionNDMap.create(
+                region=None,
+                axes=[energy],
+                unit=f_samples.unit,
+                data=np.percentile(f_samples, 84, axis=0),
+            )
+            y_lo = scale_plot_flux(y_lo, energy_power).quantity[:, 0, 0]
+            y_hi = scale_plot_flux(y_hi, energy_power).quantity[:, 0, 0]
 
         with quantity_support():
             ax.fill_between(energy.center, y_lo, y_hi, **kwargs)