(* create and initialize simulation variables ------------------------- *)
n0 = 1000                      (* number of nuclei at time t = 0        *)
tau = 1                        (* mean lifetime                         *)
tMax = 5                       (* time to stop simulation               *)
dt = 0.2                       (* time step Delta t for Euler algorithm *)
dataSets = { }                 (* list of computational results         *)
dtValues = { }                 (* list of dt values in dataSets         *)

(* define functions ----------------------------------------------------*)

getYesNo[prompt_] :=           (* function to get Yes/No response       *)
    Module[ { done = False, response, yes },
        While[ !done,
            response = InputString[prompt " Enter Y[es] or N[o]: "];
            Switch[ StringTake[ToUpperCase[response], 1],
                "Y", done = yes = True,
                "N", done = True; yes = False,
                 _ , done = False; Beep[]
            ]
        ];
        yes                    (* value returned by function            *)
    ]

getDouble[prompt_] :=          (* function to get a real number         *)
    Module[ { d },
        While[ True,
            d = Input[prompt];
            If[ NumberQ[d],
                Break[],
                Beep[]
            ]
        ];
        d                      (* value returned by function            *)
    ]

getConstants[_:Null] :=        (* user inputs N(0), tau, t_max          *)
    Module[ { },
        Print["Simulation of Radioactive Decay"];
        Print["-------------------------------"];
        Print["Current values of physical constants:"];
        Print["N(0) = ", n0, " tau = ", tau, " t_max = ", tMax];
        If[ getYesNo["Change these values?"],
            n0 = getDouble["Enter new value of N(0): "];
            tau = getDouble["Enter new value of tau: "];
            tMax = getDouble["Enter new value of t_max: "]
        ]
    ]

doComputation[_:Null] :=       (* solve equation for several dt values  *)
    Module[ { t, n, set, i },
        dataSets = { };
        dtValues = { };
        While[ True,
            dt = getDouble["Current dt = " <> ToString[dt] <>
                           ", Enter new dt:  "];
            t = 0;
            n = n0;
            set = {{t, n}};

            For[ i = 0, i < tMax/dt, i++,
                n -= n / tau * dt;
                t += dt;
                AppendTo[set, {t, n}]
            ];

            AppendTo[dataSets, set];
            AppendTo[dtValues, dt];

            If[getYesNo["Another time step?"], Continue[], Break[]]
        ]
    ]

plotResults[_:Null] :=         (* plot all results on a single graph    *)
    Module[ { t, i, dt },
        exact = Plot[n0 Exp[-t/tau], {t, 0, tMax}, DisplayFunction->Identity];
        showCommand = "Show[exact, ";
        plots = { };

        For[i = 1, i <= Length[dataSets], i++,
            plot = ListPlot[dataSets[[i]], DisplayFunction->Identity];
            AppendTo[plots, plot];
            showCommand = showCommand <> "plots[[" <> ToString[i] <> "]], "
        ];

        showCommand = showCommand <> "DisplayFunction->$DisplayFunction]";
        ToExpression[showCommand]
    ]

(* The top level program ---------------------------------------------- *)

While[ True,
    getConstants[];
    doComputation[];
    plotResults[];
    If[getYesNo["Another computation?"], Continue[], Break[]]
]