diff --git a/examples/11_plurigaussian/00_simple.py b/examples/11_plurigaussian/00_simple.py
index 264c5ea3..7e4b2dfc 100644
--- a/examples/11_plurigaussian/00_simple.py
+++ b/examples/11_plurigaussian/00_simple.py
@@ -57,9 +57,16 @@
 pgs = gs.PGS(dim, [field1, field2], L)
 
 ###############################################################################
-# Finally, we can plot the PGS, but we will also show the field `L`.
+# Finally, we can plot the PGS, but we will also show the field `L` and the two
+# original Gaussian fields.
 
-fig, axs = plt.subplots(1, 2)
+fig, axs = plt.subplots(2, 2)
 
-axs[0].imshow(L, cmap="copper")
-axs[1].imshow(pgs.P, cmap="copper")
+axs[0, 0].imshow(field1, cmap="copper")
+axs[0, 1].imshow(field2, cmap="copper")
+axs[1, 0].imshow(L, cmap="copper")
+axs[1, 1].imshow(pgs.P, cmap="copper")
+
+# For more information on Plurigaussian fields and how they naturally extend
+# truncated Gaussian fields, we can recommend the book
+# [Plurigaussian Simulations in Geosciences](https://doi.org/10.1007/978-3-642-19607-2)
diff --git a/examples/11_plurigaussian/01_pgs.py b/examples/11_plurigaussian/01_pgs.py
index e231f0e9..0b0a7f9c 100644
--- a/examples/11_plurigaussian/01_pgs.py
+++ b/examples/11_plurigaussian/01_pgs.py
@@ -30,6 +30,7 @@
 model2 = gs.Gaussian(dim=dim, var=1, len_scale=[1, 20], angles=np.pi / 4)
 srf2 = gs.SRF(model2, seed=19970221)
 field2 = srf2.structured([x, y])
+field1 += 5.0
 
 ###############################################################################
 # Internally, each field's values are mapped along an axis, which can be nicely
@@ -69,14 +70,25 @@
 # First, we calculate the indices of the L field, which we need for manually
 # plotting L together with the scatter plot. Normally they are computed internally.
 
+l = np.floor(field1.min()) - 1
+h = np.ceil(field1.max()) + 1
+m = (h + l) / 2.0
+dist = max(np.abs(h - m), np.abs(l - m))
+x_mean = field1.mean()
 x_l = np.linspace(
-    np.floor(field1[0].min()) - 1,
-    np.ceil(field1[0].max()) + 1,
+    x_mean - dist,
+    x_mean + dist,
     L.shape[0],
 )
+
+l = np.floor(field1.min()) - 1
+h = np.ceil(field1.max()) + 1
+m = (h + l) / 2.0
+dist = max(np.abs(h - m), np.abs(l - m))
+y_mean = field2.mean()
 y_l = np.linspace(
-    np.floor(field1[1].min()) - 1,
-    np.ceil(field1[1].max()) + 1,
+    y_mean - dist,
+    y_mean + dist,
     L.shape[1],
 )
 
@@ -89,10 +101,14 @@
 # And now to some plotting. Unfortunately, matplotlib likes to mess around with
 # the aspect ratios of the plots, so the left panel is a bit stretched.
 
-fig, axs = plt.subplots(1, 2)
-axs[0].scatter(field1.flatten(), field2.flatten(), s=0.1, color="C0")
-axs[0].pcolormesh(x_l, y_l, L.T, alpha=0.3)
-axs[1].imshow(pgs.P, cmap="copper")
+fig, axs = plt.subplots(2, 2)
+axs[0, 0].imshow(field1, cmap="copper")
+axs[0, 1].imshow(field2, cmap="copper")
+axs[1, 0].scatter(field1.flatten(), field2.flatten(), s=0.1, color="C0")
+axs[1, 0].pcolormesh(x_l, y_l, L.T, alpha=0.3)
+
+axs[1, 1].imshow(pgs.P, cmap="copper")
+plt.show()
 
 ###############################################################################
 # The black areas show the category 0 and the orange areas show category 1. We
diff --git a/examples/11_plurigaussian/02_spatial_relations.py b/examples/11_plurigaussian/02_spatial_relations.py
index 5e9cdfdf..320e65ca 100644
--- a/examples/11_plurigaussian/02_spatial_relations.py
+++ b/examples/11_plurigaussian/02_spatial_relations.py
@@ -87,12 +87,14 @@
 pgs = gs.PGS(dim, [field1, field2], L)
 
 ###############################################################################
-# And now the plotting of `L` and the PGS.
+# And now the plotting of the two Gaussian fields, `L`, and the PGS.
 
-fig, axs = plt.subplots(1, 2)
+fig, axs = plt.subplots(2, 2)
 
-axs[0].imshow(L, cmap="copper")
-axs[1].imshow(pgs.P, cmap="copper")
+axs[0, 0].imshow(field1, cmap="copper")
+axs[0, 1].imshow(field2, cmap="copper")
+axs[1, 0].imshow(L, cmap="copper")
+axs[1, 1].imshow(pgs.P, cmap="copper")
 
 ###############################################################################
 # We can see that the two upper light and medium brown rectangles both fill up
diff --git a/examples/11_plurigaussian/03_correlations.py b/examples/11_plurigaussian/03_correlations.py
index eeef8b4e..98008a2d 100644
--- a/examples/11_plurigaussian/03_correlations.py
+++ b/examples/11_plurigaussian/03_correlations.py
@@ -53,12 +53,14 @@
 pgs = gs.PGS(dim, [field1, field2], L)
 
 ###############################################################################
-# And now the plotting of `L` and the PGS.
+# And now the plotting of the two Gaussian fields, `L`, and the PGS.
 
-fig, axs = plt.subplots(1, 2)
+fig, axs = plt.subplots(2, 2)
 
-axs[0].imshow(L, cmap="copper")
-axs[1].imshow(pgs.P, cmap="copper")
+axs[0, 0].imshow(field1, cmap="copper")
+axs[0, 1].imshow(field2, cmap="copper")
+axs[1, 0].imshow(L, cmap="copper")
+axs[1, 1].imshow(pgs.P, cmap="copper")
 plt.show()
 
 ###############################################################################