diff --git a/demos/getting_started_1.py b/demos/getting_started_1.py
index f84d5d3b150c1a875fca1b8777366214b9dd64be..b261743b0902305950bf3f99a36cb38669e207ee 100644
--- a/demos/getting_started_1.py
+++ b/demos/getting_started_1.py
@@ -51,10 +51,15 @@ def mask_to_nan(mask, field):
 
 
 if __name__ == '__main__':
+    import sys
     np.random.seed(42)
 
     # Choose space on which the signal field is defined
-    mode = 1
+    if len(sys.argv) == 2:
+        mode = int(sys.argv[1])
+    else:
+        mode = 1
+
     if mode == 0:
         # One-dimensional regular grid
         position_space = ift.RGSpace([1024])
@@ -135,6 +140,7 @@ if __name__ == '__main__':
     # Plotting
     rg = isinstance(position_space, ift.RGSpace)
     plot = ift.Plot()
+    filename = f"getting_started_1_mode_{mode}.png"
     if rg and len(position_space.shape) == 1:
         plot.add(
             [HT(MOCK_SIGNAL), GR.adjoint(data),
@@ -142,10 +148,11 @@ if __name__ == '__main__':
             label=['Mock signal', 'Data', 'Reconstruction'],
             alpha=[1, .3, 1])
         plot.add(mask_to_nan(mask, HT(m - MOCK_SIGNAL)), title='Residuals')
-        plot.output(nx=2, ny=1, xsize=10, ysize=4, title="getting_started_1")
+        plot.output(nx=2, ny=1, xsize=10, ysize=4, name=filename)
     else:
         plot.add(HT(MOCK_SIGNAL), title='Mock Signal')
         plot.add(mask_to_nan(mask, (GR(Mask)).adjoint(data)), title='Data')
         plot.add(HT(m), title='Reconstruction')
         plot.add(mask_to_nan(mask, HT(m - MOCK_SIGNAL)), title='Residuals')
-        plot.output(nx=2, ny=2, xsize=10, ysize=10, title="getting_started_1")
+        plot.output(nx=2, ny=2, xsize=10, ysize=10, name=filename)
+    print(f"Saved results as '{filename}'.")
diff --git a/demos/getting_started_2.py b/demos/getting_started_2.py
index 5374731bc0629c8d6f3f62fdc2491d6fae6e28b4..a1acc7296e2cab96502b6bc0c6e4034abc22b3be 100644
--- a/demos/getting_started_2.py
+++ b/demos/getting_started_2.py
@@ -42,11 +42,16 @@ def exposure_2d():
 
 
 if __name__ == '__main__':
+    import sys
     # FIXME All random seeds to 42
     np.random.seed(42)
 
     # Choose space on which the signal field is defined
-    mode = 1
+    if len(sys.argv) == 2:
+        mode = int(sys.argv[1])
+    else:
+        mode = 1
+
     if mode == 0:
         # One-dimensional regular grid with uniform exposure of 10
         position_space = ift.RGSpace(1024)
@@ -107,9 +112,11 @@ if __name__ == '__main__':
     # Plotting
     signal = sky(mock_position)
     reconst = sky(H.position)
+    filename = f"getting_started_2_mode_{mode}.png"
     plot = ift.Plot()
     plot.add(signal, title='Signal')
     plot.add(GR.adjoint(data), title='Data')
     plot.add(reconst, title='Reconstruction')
     plot.add(reconst - signal, title='Residuals')
-    plot.output(name='getting_started_2.pdf', xsize=16, ysize=16)
+    plot.output(xsize=12, ysize=10, name=filename)
+    print(f"Saved results as '{filename}'.")
diff --git a/demos/getting_started_3.py b/demos/getting_started_3.py
index 59145c408103442ef586dce07997b7b0eeb3add1..c69f059c5f55e94570dc66ef9f8bbcfc9d8b6e8b 100644
--- a/demos/getting_started_3.py
+++ b/demos/getting_started_3.py
@@ -43,11 +43,16 @@ def radial_los(n_los):
 
 
 if __name__ == '__main__':
+    import sys
     np.random.seed(420)  # picked for a nice field realization
 
     # Choose between random line-of-sight response (mode=0) and radial lines
     # of sight (mode=1)
-    mode = 0
+    if len(sys.argv) == 2:
+        mode = int(sys.argv[1])
+    else:
+        mode = 0
+    filename = f"getting_started_3_mode_{mode}_" + "{}.png"
 
     position_space = ift.RGSpace([128, 128])
     harmonic_space = position_space.get_default_codomain()
@@ -113,7 +118,7 @@ if __name__ == '__main__':
     plot.add(signal(mock_position), title='Ground Truth')
     plot.add(R.adjoint_times(data), title='Data')
     plot.add([A.force(mock_position)], title='Power Spectrum')
-    plot.output(ny=1, nx=3, xsize=24, ysize=6, name="setup.png")
+    plot.output(ny=1, nx=3, xsize=24, ysize=6, name=filename.format("setup"))
 
     # number of samples used to estimate the KL
     N_samples = 20
@@ -129,7 +134,8 @@ if __name__ == '__main__':
         plot = ift.Plot()
         plot.add(signal(KL.position), title="reconstruction")
         plot.add([A.force(KL.position), A.force(mock_position)], title="power")
-        plot.output(ny=1, ysize=6, xsize=16, name="loop-{:02}.png".format(i))
+        plot.output(ny=1, ysize=6, xsize=16,
+                    name=filename.format(f"loop_{i:02}"))
 
     # Draw posterior samples
     KL = ift.MetricGaussianKL(mean, H, N_samples)
@@ -138,6 +144,7 @@ if __name__ == '__main__':
         sc.add(signal(sample + KL.position))
 
     # Plotting
+    filename_res = filename.format("results")
     plot = ift.Plot()
     plot.add(sc.mean, title="Posterior Mean")
     plot.add(ift.sqrt(sc.var), title="Posterior Standard Deviation")
@@ -148,4 +155,5 @@ if __name__ == '__main__':
                   A.force(mock_position)],
         title="Sampled Posterior Power Spectrum",
         linewidth=[1.]*len(powers) + [3., 3.])
-    plot.output(ny=1, nx=3, xsize=24, ysize=6, name="results.png")
+    plot.output(ny=1, nx=3, xsize=24, ysize=6, name=filename_res)
+    print(f"Saved results as '{filename_res}'.")