diff --git a/src/core/examples/main-random-variable-stream.cc b/src/core/examples/main-random-variable-stream.cc
index b0455524d..7caea9925 100644
--- a/src/core/examples/main-random-variable-stream.cc
+++ b/src/core/examples/main-random-variable-stream.cc
@@ -25,6 +25,8 @@
 #include "ns3/string.h"
 
 #include <cmath>
+#include <fstream>
+#include <iostream>
 #include <map>
 
 /**
@@ -116,15 +118,17 @@ main(int argc, char* argv[])
     unsigned int probes = 1000000;
     double precision = 0.01;
 
-    GnuplotCollection gnuplots("main-random-variables.pdf");
+    const std::string plotFile{"main-random-variables"};
+    GnuplotCollection gnuplots(plotFile + ".pdf");
     gnuplots.SetTerminal("pdf enhanced");
 
     {
+        std::cout << "UniformRandomVariable......." << std::flush;
         Gnuplot plot;
         plot.SetTitle("UniformRandomVariable");
         plot.AppendExtra("set yrange [0:]");
 
-        Ptr<UniformRandomVariable> x = CreateObject<UniformRandomVariable>();
+        auto x = CreateObject<UniformRandomVariable>();
         x->SetAttribute("Min", DoubleValue(0.0));
         x->SetAttribute("Max", DoubleValue(1.0));
 
@@ -132,29 +136,31 @@ main(int argc, char* argv[])
         plot.AddDataset(Gnuplot2dFunction("1.0", "0 <= x && x <= 1 ? 1.0 : 0"));
 
         gnuplots.AddPlot(plot);
+        std::cout << "done" << std::endl;
     }
 
     {
+        std::cout << "ExponentialRandomVariable..." << std::flush;
         Gnuplot plot;
         plot.SetTitle("ExponentialRandomVariable");
-        plot.AppendExtra("set xrange [0:8]");
+        plot.AppendExtra("set xrange [0:4]");
         plot.AppendExtra("ExpDist(x,l) = 1/l * exp(-1/l * x)");
 
-        Ptr<ExponentialRandomVariable> x1 = CreateObject<ExponentialRandomVariable>();
+        auto x1 = CreateObject<ExponentialRandomVariable>();
         x1->SetAttribute("Mean", DoubleValue(0.5));
 
         plot.AddDataset(Histogram(x1, probes, precision, "ExponentialRandomVariable m=0.5"));
 
         plot.AddDataset(Gnuplot2dFunction("ExponentialDistribution mean 0.5", "ExpDist(x, 0.5)"));
 
-        Ptr<ExponentialRandomVariable> x2 = CreateObject<ExponentialRandomVariable>();
+        auto x2 = CreateObject<ExponentialRandomVariable>();
         x2->SetAttribute("Mean", DoubleValue(1.0));
 
         plot.AddDataset(Histogram(x2, probes, precision, "ExponentialRandomVariable m=1"));
 
         plot.AddDataset(Gnuplot2dFunction("ExponentialDistribution mean 1.0", "ExpDist(x, 1.0)"));
 
-        Ptr<ExponentialRandomVariable> x3 = CreateObject<ExponentialRandomVariable>();
+        auto x3 = CreateObject<ExponentialRandomVariable>();
         x3->SetAttribute("Mean", DoubleValue(1.5));
 
         plot.AddDataset(Histogram(x3, probes, precision, "ExponentialRandomVariable m=1.5"));
@@ -162,28 +168,30 @@ main(int argc, char* argv[])
         plot.AddDataset(Gnuplot2dFunction("ExponentialDistribution mean 1.5", "ExpDist(x, 1.5)"));
 
         gnuplots.AddPlot(plot);
+        std::cout << "done" << std::endl;
     }
 
     {
+        std::cout << "ParetoRandomVariable........" << std::flush;
         Gnuplot plot;
         plot.SetTitle("ParetoRandomVariable");
-        plot.AppendExtra("set xrange [0:2]");
+        plot.AppendExtra("set xrange [0:4]");
 
-        Ptr<ParetoRandomVariable> x1 = CreateObject<ParetoRandomVariable>();
+        auto x1 = CreateObject<ParetoRandomVariable>();
         x1->SetAttribute("Scale", DoubleValue(1.0));
         x1->SetAttribute("Shape", DoubleValue(1.5));
 
         plot.AddDataset(
             Histogram(x1, probes, precision, "ParetoRandomVariable scale=1.0 shape=1.5"));
 
-        Ptr<ParetoRandomVariable> x2 = CreateObject<ParetoRandomVariable>();
+        auto x2 = CreateObject<ParetoRandomVariable>();
         x2->SetAttribute("Scale", DoubleValue(1.0));
         x2->SetAttribute("Shape", DoubleValue(2.0));
 
         plot.AddDataset(
             Histogram(x2, probes, precision, "ParetoRandomVariable scale=1.0 shape=2.0"));
 
-        Ptr<ParetoRandomVariable> x3 = CreateObject<ParetoRandomVariable>();
+        auto x3 = CreateObject<ParetoRandomVariable>();
         x3->SetAttribute("Scale", DoubleValue(1.0));
         x3->SetAttribute("Shape", DoubleValue(2.5));
 
@@ -191,28 +199,30 @@ main(int argc, char* argv[])
             Histogram(x3, probes, precision, "ParetoRandomVariable scale=1.0 shape=2.5"));
 
         gnuplots.AddPlot(plot);
+        std::cout << "done" << std::endl;
     }
 
     {
+        std::cout << "WeibullRandomVariable......." << std::flush;
         Gnuplot plot;
         plot.SetTitle("WeibullRandomVariable");
-        plot.AppendExtra("set xrange [0:3]");
+        plot.AppendExtra("set xrange [0:4]");
 
-        Ptr<WeibullRandomVariable> x1 = CreateObject<WeibullRandomVariable>();
+        auto x1 = CreateObject<WeibullRandomVariable>();
         x1->SetAttribute("Scale", DoubleValue(1.0));
         x1->SetAttribute("Shape", DoubleValue(1.0));
 
         plot.AddDataset(
             Histogram(x1, probes, precision, "WeibullRandomVariable scale=1.0 shape=1.0"));
 
-        Ptr<WeibullRandomVariable> x2 = CreateObject<WeibullRandomVariable>();
+        auto x2 = CreateObject<WeibullRandomVariable>();
         x2->SetAttribute("Scale", DoubleValue(1.0));
         x2->SetAttribute("Shape", DoubleValue(2.0));
 
         plot.AddDataset(
             Histogram(x2, probes, precision, "WeibullRandomVariable scale=1.0 shape=2.0"));
 
-        Ptr<WeibullRandomVariable> x3 = CreateObject<WeibullRandomVariable>();
+        auto x3 = CreateObject<WeibullRandomVariable>();
         x3->SetAttribute("Scale", DoubleValue(1.0));
         x3->SetAttribute("Shape", DoubleValue(3.0));
 
@@ -220,16 +230,18 @@ main(int argc, char* argv[])
             Histogram(x3, probes, precision, "WeibullRandomVariable scale=1.0 shape=3.0"));
 
         gnuplots.AddPlot(plot);
+        std::cout << "done" << std::endl;
     }
 
     {
+        std::cout << "NormalRandomVariable........" << std::flush;
         Gnuplot plot;
         plot.SetTitle("NormalRandomVariable");
-        plot.AppendExtra("set xrange [-3:3]");
+        plot.AppendExtra("set xrange [-4:4]");
         plot.AppendExtra(
             "NormalDist(x,m,s) = 1 / (s * sqrt(2*pi)) * exp(-1.0 / 2.0 * ((x-m) / s)**2)");
 
-        Ptr<NormalRandomVariable> x1 = CreateObject<NormalRandomVariable>();
+        auto x1 = CreateObject<NormalRandomVariable>();
         x1->SetAttribute("Mean", DoubleValue(0.0));
         x1->SetAttribute("Variance", DoubleValue(1.0));
 
@@ -238,7 +250,7 @@ main(int argc, char* argv[])
         plot.AddDataset(Gnuplot2dFunction("NormalDist {/Symbol m}=0.0 {/Symbol s}=1.0",
                                           "NormalDist(x,0.0,1.0)"));
 
-        Ptr<NormalRandomVariable> x2 = CreateObject<NormalRandomVariable>();
+        auto x2 = CreateObject<NormalRandomVariable>();
         x2->SetAttribute("Mean", DoubleValue(0.0));
         x2->SetAttribute("Variance", DoubleValue(2.0));
 
@@ -247,7 +259,7 @@ main(int argc, char* argv[])
         plot.AddDataset(Gnuplot2dFunction("NormalDist {/Symbol m}=0.0 {/Symbol s}=sqrt(2.0)",
                                           "NormalDist(x,0.0,sqrt(2.0))"));
 
-        Ptr<NormalRandomVariable> x3 = CreateObject<NormalRandomVariable>();
+        auto x3 = CreateObject<NormalRandomVariable>();
         x3->SetAttribute("Mean", DoubleValue(0.0));
         x3->SetAttribute("Variance", DoubleValue(3.0));
 
@@ -257,55 +269,61 @@ main(int argc, char* argv[])
                                           "NormalDist(x,0.0,sqrt(3.0))"));
 
         gnuplots.AddPlot(plot);
+        std::cout << "done" << std::endl;
     }
 
-    /** \todo Turn this plot back on once its distribution has been finished. */
-    /*
     {
-      Gnuplot plot;
-      plot.SetTitle ("EmpiricalRandomVariable");
-      plot.AppendExtra ("set xrange [*:*]");
+        std::cout << "EmpiricalVariable..........." << std::flush;
+        Gnuplot plot;
+        plot.SetTitle("EmpiricalRandomVariable");
+        plot.AppendExtra("set xrange [*:*]");
 
-      EmpiricalRandomVariable emp1;
-      emp1.CDF (0.0,  0.0 / 15.0);
-      emp1.CDF (0.2,  1.0 / 15.0);
-      emp1.CDF (0.4,  3.0 / 15.0);
-      emp1.CDF (0.6,  6.0 / 15.0);
-      emp1.CDF (0.8, 10.0 / 15.0);
-      emp1.CDF (1.0, 15.0 / 15.0);
+        auto x = CreateObject<EmpiricalRandomVariable>();
+        x->CDF(0.0, 0.0 / 15.0);
+        x->CDF(0.2, 1.0 / 15.0);
+        x->CDF(0.4, 3.0 / 15.0);
+        x->CDF(0.6, 6.0 / 15.0);
+        x->CDF(0.8, 10.0 / 15.0);
+        x->CDF(1.0, 15.0 / 15.0);
 
-      plot.AddDataset ( Histogram (emp1, probes, precision,
-                                    "EmpiricalRandomVariable (Stairs)") );
+        plot.AddDataset(
+            Histogram(x, probes, precision, "EmpiricalRandomVariable (Sampling)", true));
 
-      gnuplots.AddPlot (plot);
+        x->SetInterpolate(true);
+        plot.AppendExtra("set y2range [0:*]");
+        auto d2 = Histogram(x, probes, precision, "EmpiricalRandomVariable (Interpolate)");
+        d2.SetExtra(" axis x1y2");
+        plot.AddDataset(d2);
+
+        gnuplots.AddPlot(plot);
+        std::cout << "done" << std::endl;
     }
-    */
 
-    /** \todo Turn this plot back on once its distribution has been finished. */
-    /*
     {
-      Gnuplot plot;
-      plot.SetTitle ("DeterministicRandomVariable");
-      plot.AppendExtra ("set xrange [*:*]");
+        std::cout << "DeterministicVariable......." << std::flush;
+        Gnuplot plot;
+        plot.SetTitle("DeterministicRandomVariable");
+        plot.AppendExtra("set xrange [*:*]");
 
-      double values[] = { 0.0, 0.2, 0.2, 0.4, 0.2, 0.6, 0.8, 0.8, 1.0 };
-      DeterministicRandomVariable det1 (values, sizeof(values) / sizeof(values[0]));
+        auto x1 = CreateObject<DeterministicRandomVariable>();
+        double values[] = {0.0, 0.2, 0.2, 0.4, 0.2, 0.6, 0.8, 0.8, 1.0};
+        x1->SetValueArray(values, sizeof(values) / sizeof(values[0]));
 
-      plot.AddDataset ( Histogram (det1, probes, precision,
-                                    "DeterministicRandomVariable", true) );
+        plot.AddDataset(Histogram(x1, probes, precision, "DeterministicRandomVariable", true));
 
-      gnuplots.AddPlot (plot);
+        gnuplots.AddPlot(plot);
+        std::cout << "done" << std::endl;
     }
-    */
 
     {
+        std::cout << "LogNormalRandomVariable....." << std::flush;
         Gnuplot plot;
         plot.SetTitle("LogNormalRandomVariable");
-        plot.AppendExtra("set xrange [0:3]");
+        plot.AppendExtra("set xrange [0:4]");
 
         plot.AppendExtra("LogNormalDist(x,m,s) = 1.0/x * NormalDist(log(x), m, s)");
 
-        Ptr<LogNormalRandomVariable> x1 = CreateObject<LogNormalRandomVariable>();
+        auto x1 = CreateObject<LogNormalRandomVariable>();
         x1->SetAttribute("Mu", DoubleValue(0.0));
         x1->SetAttribute("Sigma", DoubleValue(1.0));
 
@@ -314,13 +332,13 @@ main(int argc, char* argv[])
         plot.AddDataset(
             Gnuplot2dFunction("LogNormalDist(x, 0.0, 1.0)", "LogNormalDist(x, 0.0, 1.0)"));
 
-        Ptr<LogNormalRandomVariable> x2 = CreateObject<LogNormalRandomVariable>();
+        auto x2 = CreateObject<LogNormalRandomVariable>();
         x2->SetAttribute("Mu", DoubleValue(0.0));
         x2->SetAttribute("Sigma", DoubleValue(0.5));
 
         plot.AddDataset(Histogram(x2, probes, precision, "LogNormalRandomVariable m=0.0 s=0.5"));
 
-        Ptr<LogNormalRandomVariable> x3 = CreateObject<LogNormalRandomVariable>();
+        auto x3 = CreateObject<LogNormalRandomVariable>();
         x3->SetAttribute("Mu", DoubleValue(0.0));
         x3->SetAttribute("Sigma", DoubleValue(0.25));
 
@@ -329,13 +347,13 @@ main(int argc, char* argv[])
         plot.AddDataset(
             Gnuplot2dFunction("LogNormalDist(x, 0.0, 0.25)", "LogNormalDist(x, 0.0, 0.25)"));
 
-        Ptr<LogNormalRandomVariable> x4 = CreateObject<LogNormalRandomVariable>();
+        auto x4 = CreateObject<LogNormalRandomVariable>();
         x4->SetAttribute("Mu", DoubleValue(0.0));
         x4->SetAttribute("Sigma", DoubleValue(0.125));
 
         plot.AddDataset(Histogram(x4, probes, precision, "LogNormalRandomVariable m=0.0 s=0.125"));
 
-        Ptr<LogNormalRandomVariable> x5 = CreateObject<LogNormalRandomVariable>();
+        auto x5 = CreateObject<LogNormalRandomVariable>();
         x5->SetAttribute("Mu", DoubleValue(0.0));
         x5->SetAttribute("Sigma", DoubleValue(2.0));
 
@@ -344,21 +362,23 @@ main(int argc, char* argv[])
         plot.AddDataset(
             Gnuplot2dFunction("LogNormalDist(x, 0.0, 2.0)", "LogNormalDist(x, 0.0, 2.0)"));
 
-        Ptr<LogNormalRandomVariable> x6 = CreateObject<LogNormalRandomVariable>();
+        auto x6 = CreateObject<LogNormalRandomVariable>();
         x6->SetAttribute("Mu", DoubleValue(0.0));
         x6->SetAttribute("Sigma", DoubleValue(2.5));
 
         plot.AddDataset(Histogram(x6, probes, precision, "LogNormalRandomVariable m=0.0 s=2.5"));
 
         gnuplots.AddPlot(plot);
+        std::cout << "done" << std::endl;
     }
 
     {
+        std::cout << "TriangularRandomVariable...." << std::flush;
         Gnuplot plot;
         plot.SetTitle("TriangularRandomVariable");
         plot.AppendExtra("set xrange [*:*]");
 
-        Ptr<TriangularRandomVariable> x1 = CreateObject<TriangularRandomVariable>();
+        auto x1 = CreateObject<TriangularRandomVariable>();
         x1->SetAttribute("Min", DoubleValue(0.0));
         x1->SetAttribute("Max", DoubleValue(1.0));
         x1->SetAttribute("Mean", DoubleValue(0.5));
@@ -366,7 +386,7 @@ main(int argc, char* argv[])
         plot.AddDataset(
             Histogram(x1, probes, precision, "TriangularRandomVariable [0.0 .. 1.0) m=0.5"));
 
-        Ptr<TriangularRandomVariable> x2 = CreateObject<TriangularRandomVariable>();
+        auto x2 = CreateObject<TriangularRandomVariable>();
         x2->SetAttribute("Min", DoubleValue(0.0));
         x2->SetAttribute("Max", DoubleValue(1.0));
         x2->SetAttribute("Mean", DoubleValue(0.4));
@@ -374,7 +394,7 @@ main(int argc, char* argv[])
         plot.AddDataset(
             Histogram(x2, probes, precision, "TriangularRandomVariable [0.0 .. 1.0) m=0.4"));
 
-        Ptr<TriangularRandomVariable> x3 = CreateObject<TriangularRandomVariable>();
+        auto x3 = CreateObject<TriangularRandomVariable>();
         x3->SetAttribute("Min", DoubleValue(0.0));
         x3->SetAttribute("Max", DoubleValue(1.0));
         x3->SetAttribute("Mean", DoubleValue(0.65));
@@ -383,9 +403,11 @@ main(int argc, char* argv[])
             Histogram(x3, probes, precision, "TriangularRandomVariable [0.0 .. 1.0) m=0.65"));
 
         gnuplots.AddPlot(plot);
+        std::cout << "done" << std::endl;
     }
 
     {
+        std::cout << "GammaRandomVariable........." << std::flush;
         Gnuplot plot;
         plot.SetTitle("GammaRandomVariable");
         plot.AppendExtra("set xrange [0:10]");
@@ -396,7 +418,7 @@ main(int argc, char* argv[])
             "set label 1 '{/Symbol g}(x,{/Symbol a},{/Symbol b}) = x^{/Symbol a-1} e^{-x {/Symbol "
             "b}^{-1}} ( {/Symbol b}^{/Symbol a} {/Symbol G}({/Symbol a}) )^{-1}' at 0.7, 0.9");
 
-        Ptr<GammaRandomVariable> x1 = CreateObject<GammaRandomVariable>();
+        auto x1 = CreateObject<GammaRandomVariable>();
         x1->SetAttribute("Alpha", DoubleValue(1.0));
         x1->SetAttribute("Beta", DoubleValue(1.0));
 
@@ -404,7 +426,7 @@ main(int argc, char* argv[])
 
         plot.AddDataset(Gnuplot2dFunction("{/Symbol g}(x, 1.0, 1.0)", "GammaDist(x, 1.0, 1.0)"));
 
-        Ptr<GammaRandomVariable> x2 = CreateObject<GammaRandomVariable>();
+        auto x2 = CreateObject<GammaRandomVariable>();
         x2->SetAttribute("Alpha", DoubleValue(1.5));
         x2->SetAttribute("Beta", DoubleValue(1.0));
 
@@ -412,7 +434,7 @@ main(int argc, char* argv[])
 
         plot.AddDataset(Gnuplot2dFunction("{/Symbol g}(x, 1.5, 1.0)", "GammaDist(x, 1.5, 1.0)"));
 
-        Ptr<GammaRandomVariable> x3 = CreateObject<GammaRandomVariable>();
+        auto x3 = CreateObject<GammaRandomVariable>();
         x3->SetAttribute("Alpha", DoubleValue(2.0));
         x3->SetAttribute("Beta", DoubleValue(1.0));
 
@@ -420,7 +442,7 @@ main(int argc, char* argv[])
 
         plot.AddDataset(Gnuplot2dFunction("{/Symbol g}(x, 2.0, 1.0)", "GammaDist(x, 2.0, 1.0)"));
 
-        Ptr<GammaRandomVariable> x4 = CreateObject<GammaRandomVariable>();
+        auto x4 = CreateObject<GammaRandomVariable>();
         x4->SetAttribute("Alpha", DoubleValue(4.0));
         x4->SetAttribute("Beta", DoubleValue(1.0));
 
@@ -428,7 +450,7 @@ main(int argc, char* argv[])
 
         plot.AddDataset(Gnuplot2dFunction("{/Symbol g}(x, 4.0, 1.0)", "GammaDist(x, 4.0, 1.0)"));
 
-        Ptr<GammaRandomVariable> x5 = CreateObject<GammaRandomVariable>();
+        auto x5 = CreateObject<GammaRandomVariable>();
         x5->SetAttribute("Alpha", DoubleValue(2.0));
         x5->SetAttribute("Beta", DoubleValue(2.0));
 
@@ -436,7 +458,7 @@ main(int argc, char* argv[])
 
         plot.AddDataset(Gnuplot2dFunction("{/Symbol g}(x, 2.0, 2.0)", "GammaDist(x, 2.0, 2.0)"));
 
-        Ptr<GammaRandomVariable> x6 = CreateObject<GammaRandomVariable>();
+        auto x6 = CreateObject<GammaRandomVariable>();
         x6->SetAttribute("Alpha", DoubleValue(2.5));
         x6->SetAttribute("Beta", DoubleValue(3.0));
 
@@ -444,7 +466,7 @@ main(int argc, char* argv[])
 
         plot.AddDataset(Gnuplot2dFunction("{/Symbol g}(x, 2.5, 3.0)", "GammaDist(x, 2.5, 3.0)"));
 
-        Ptr<GammaRandomVariable> x7 = CreateObject<GammaRandomVariable>();
+        auto x7 = CreateObject<GammaRandomVariable>();
         x7->SetAttribute("Alpha", DoubleValue(2.5));
         x7->SetAttribute("Beta", DoubleValue(4.5));
 
@@ -453,9 +475,11 @@ main(int argc, char* argv[])
         plot.AddDataset(Gnuplot2dFunction("{/Symbol g}(x, 2.5, 4.5)", "GammaDist(x, 2.5, 4.5)"));
 
         gnuplots.AddPlot(plot);
+        std::cout << "done" << std::endl;
     }
 
     {
+        std::cout << "ErlangRandomVariable........" << std::flush;
         Gnuplot plot;
         plot.SetTitle("ErlangRandomVariable");
         plot.AppendExtra("set xrange [0:10]");
@@ -464,7 +488,7 @@ main(int argc, char* argv[])
         plot.AppendExtra("set label 1 'Erlang(x,k,{/Symbol l}) = x^{k-1} e^{-x {/Symbol l}^{-1}} ( "
                          "{/Symbol l}^k (k-1)! )^{-1}' at 0.7, 0.9");
 
-        Ptr<ErlangRandomVariable> x1 = CreateObject<ErlangRandomVariable>();
+        auto x1 = CreateObject<ErlangRandomVariable>();
         x1->SetAttribute("K", IntegerValue(1));
         x1->SetAttribute("Lambda", DoubleValue(1.0));
 
@@ -473,7 +497,7 @@ main(int argc, char* argv[])
 
         plot.AddDataset(Gnuplot2dFunction("Erlang(x, 1, 1.0)", "ErlangDist(x, 1, 1.0)"));
 
-        Ptr<ErlangRandomVariable> x2 = CreateObject<ErlangRandomVariable>();
+        auto x2 = CreateObject<ErlangRandomVariable>();
         x2->SetAttribute("K", IntegerValue(2));
         x2->SetAttribute("Lambda", DoubleValue(1.0));
 
@@ -482,7 +506,7 @@ main(int argc, char* argv[])
 
         plot.AddDataset(Gnuplot2dFunction("Erlang(x, 2, 1.0)", "ErlangDist(x, 2, 1.0)"));
 
-        Ptr<ErlangRandomVariable> x3 = CreateObject<ErlangRandomVariable>();
+        auto x3 = CreateObject<ErlangRandomVariable>();
         x3->SetAttribute("K", IntegerValue(3));
         x3->SetAttribute("Lambda", DoubleValue(1.0));
 
@@ -491,7 +515,7 @@ main(int argc, char* argv[])
 
         plot.AddDataset(Gnuplot2dFunction("Erlang(x, 3, 1.0)", "ErlangDist(x, 3, 1.0)"));
 
-        Ptr<ErlangRandomVariable> x4 = CreateObject<ErlangRandomVariable>();
+        auto x4 = CreateObject<ErlangRandomVariable>();
         x4->SetAttribute("K", IntegerValue(5));
         x4->SetAttribute("Lambda", DoubleValue(1.0));
 
@@ -500,7 +524,7 @@ main(int argc, char* argv[])
 
         plot.AddDataset(Gnuplot2dFunction("Erlang(x, 5, 1.0)", "ErlangDist(x, 5, 1.0)"));
 
-        Ptr<ErlangRandomVariable> x5 = CreateObject<ErlangRandomVariable>();
+        auto x5 = CreateObject<ErlangRandomVariable>();
         x5->SetAttribute("K", IntegerValue(2));
         x5->SetAttribute("Lambda", DoubleValue(2.0));
 
@@ -509,7 +533,7 @@ main(int argc, char* argv[])
 
         plot.AddDataset(Gnuplot2dFunction("Erlang(x, 2, 2.0)", "ErlangDist(x, 2, 2.0)"));
 
-        Ptr<ErlangRandomVariable> x6 = CreateObject<ErlangRandomVariable>();
+        auto x6 = CreateObject<ErlangRandomVariable>();
         x6->SetAttribute("K", IntegerValue(2));
         x6->SetAttribute("Lambda", DoubleValue(3.0));
 
@@ -518,7 +542,7 @@ main(int argc, char* argv[])
 
         plot.AddDataset(Gnuplot2dFunction("Erlang(x, 2, 3.0)", "ErlangDist(x, 2, 3.0)"));
 
-        Ptr<ErlangRandomVariable> x7 = CreateObject<ErlangRandomVariable>();
+        auto x7 = CreateObject<ErlangRandomVariable>();
         x7->SetAttribute("K", IntegerValue(2));
         x7->SetAttribute("Lambda", DoubleValue(5.0));
 
@@ -528,9 +552,21 @@ main(int argc, char* argv[])
         plot.AddDataset(Gnuplot2dFunction("Erlang(x, 2, 5.0)", "ErlangDist(x, 2, 5.0)"));
 
         gnuplots.AddPlot(plot);
+        std::cout << "done" << std::endl;
     }
 
-    gnuplots.GenerateOutput(std::cout);
+    {
+        std::string gnuFile = plotFile + ".plt";
+        std::cout << "Writing Gnuplot file: " << gnuFile << "..." << std::flush;
+        std::ofstream gnuStream(gnuFile);
+        gnuplots.GenerateOutput(gnuStream);
+        gnuStream.close();
+
+        std::cout << "done\nGenerating " << plotFile << ".pdf..." << std::flush;
+        std::string shellcmd = "gnuplot " + gnuFile;
+        std::system(shellcmd.c_str());
+        std::cout << "done" << std::endl;
+    }
 
     return 0;
 }