How to determine how many runs are needed for a Monte Carlo simulation? I have been reading a lot of references online, but people seem to have different sorts of equations to determine the required number of runs.

This formula was the one I obtained from an article and I am wondering if this is indeed the right one. Also, how can I use this formula?
Edit:
I will use the Monte Carlo simulation to generate random numbers from a normal distribution. These random numbers will be used to represent demand as the random variates.
The model of interest is an inventory model where demand will be the variable to be simulated. The properties of the demand is that it is normally distributed, mean is non-stationary over time, independent and not correlated, and stochastic over a time horizon.
I performed a linear program using MATLAB, however, the solution of the problem is deterministic by nature--meaning it will always give the same "optimum" value for the specified parameter settings. What I want to achieve is to incorporate randomness into the system by using a Monte Carlo simulation which will test my model for robustness when different "real" demands are incorporated.
I would say a 95-98% accuracy would be good, but it is completely arbitrary because I do not really have anything to base it on. 
Below is my code for the generation of the "real" demand values.
Where:
n= number of runs;
average(i)=mean of demand;
a(i)=standard deviation of demand;
CV=coefficient of variation;
sim_d="real" demand values;

This is the flowchart of the process to further explain the procedure (click to see)

 A: The setting is a bit different from what I'm used to, but I think this is a classic Uncertainty Quantification (UQ) approach: specifically, it's Forward Uncertainty Propagation. The approach has limits, since you assume the demand distribution, instead than estimating it from data. Anyway, let's consider the following setting: with your model, you compute $C$, the number of burgers ordered over $p$ days. Using a constant coefficient of variation $\lambda$ and a vector of means $\boldsymbol{\mu}=(\mu_1,\dots, \mu_p)$, you define a random vector $\mathbf{d}=(d_1,\dots,d_p)$, whose components are independent, normally distributed with means given by $\boldsymbol{\mu}$ and standard deviations $(\sigma_1=\lambda\mu_1,\dots, \sigma_p=\lambda\mu_p)$. Then, we generate a random sample $D=(\mathbf{d}_1,\dots,\mathbf{d}_N)$ of size $N$ from the distribution of this vector:
set.seed(25)
C <- 50
N <- 500
lambda <- 0.1
means <- c(30, 20)
p <- length(means)

demands <- rnorm(N * p, mean = means, sd = lambda * means ) # vector of length p*N

Now, for each vector $\mathbf{d}_i=(d_{i1},\dots,d_{ip})$, you either had a stock-out event or not, depending on whether $\sum_{j=1}^p d_{ij} > C$ or not. In other words, the occurence of a stock-out is a Bernoulli random variable with mean $p$ and standard deviation $(1-p)p$. Let's build the vector of stock-out events:
seqs <- seq_along(demands)
orders_over_p_days  <- tapply(demands, rep(seqs, each = p)[seqs], FUN = sum) # vector of length N
stock_outs <- orders_over_p_days > C # logical vector, 0 if we had enough burgers for p days otherwise 1

The Monte Carlo method estimates the probability $p$ of a stock-out event using the sample mean of this vector. In other words:
$$\hat{p}_N=\frac{\sum_{i=1}^N p_i}{N}$$
$\hat{p}_N$ is the Monte Carlo estimator. Note that for each $i$, $p_i$ is either 0 (no stock-out) or 1 (stock-out). For example, in our case the value of the Monte carlo estimator is
phat <- sum(stock_outs)/N
> phat
#[1] 0.488

Since we have independent and identically distributed samples $p_i$, then, for $N\to\infty$, the distribution of this random variable
$$Z_N=\sqrt{N}\frac{\hat{p}_N-p}{\sqrt{(1-p)p}}$$
converges pointwise to a standard normal distribution. You can use this result to get approximate confidence intervals for your estimate of $p$, and from that, an estimate for the number of runs required to estimate $p$ to within a tolerance $\epsilon$, with a confidence level $\gamma$. Skipping the math, an approximate expression for $N$ would be
$$N=\frac{(1-p)p\left(\Phi^{-1}(\frac{1+\gamma}{2})\right)^2}{\epsilon^2}$$
where $\Phi^{-1}(x)$ is the quantile function of the standard normal distribution. Clearly, there is a problem here: the estimate for $N$ contains the unknown $p$. A practical solution is to run the Monte Carlo simulation for an initial number of runs, say, $N_{in}=500$, and compute $\hat{p}_{N_{in}}$. Using $\hat{p}_{N_{in}}$ we get an initial, crude estimate $(1-\hat{p}_{N_{in}})\hat{p}_{N_{in}}$ for $(1-p)p$, which leads to the approximate estimate for the total number of runs:
$$N=\frac{(1-\hat{p}_{N_{in}})\hat{p}_{N_{in}}\left(\Phi^{-1}(\frac{1+\gamma}{2})\right)^2}{\epsilon^2}$$
In deriving this formula we implicitly used quite a crude estimate for the confidence interval. For this reason, it's a good idea to increase the number of runs with respect to what you got from this formula, and to check the behavior of your estimator as the number of runs increases.
As an illustration, let's assume a confidence level $\gamma=0.95$ and an absolute accuracy $\epsilon = 0.01$, meaning that you would like to estimate the probability of a stock-out in $p$ days with an error of $\pm 0.01$, at a confidence level $0.95$. In our example we get:
gamma   <- 0.95
epsilon <- 0.01
q <- qnorm((1+gamma)/2)
N <- ceiling((1-phat)*phat*q^2/epsilon^2)
#[1] 9599

I would usually increase $N$ consistently, but let's use $N=10000$ just for the sake of example:
N <- 10000
demands <- rnorm(N * p, mean = means, sd = lambda * means )    
seqs <- seq_along(demands)
orders_over_p_days  <- tapply(demands, rep(seqs, each = p)[seqs], FUN = sum)     
stock_outs <- C < orders_over_p_days 

phat <- sum(stock_outs)/N
#[1] 0.4923

We can also plot the trend of our estimate against $N$:
sample_size <- seq_along(stock_outs)
phat <- cumsum(stock_outs)/sample_size
op <- par(mar = c(5, 5, 4, 2) + 0.1)
plot(sample_size, phat, type = "l", xlab = "N", ylab =expression(hat(p)))


We conclude that an (approximate) $95\%-$confidence interval for $p$ is $0.49\pm0.01$.
PS the above formula is valid if $\epsilon$ is an absolute error: if you're interested in the relative error, then the approximate formula becomes
$$N=\frac{\frac{(1-\hat{p}_{N_{in}})}{\hat{p}_{N_{in}}}\left(\Phi^{-1}(\frac{1+\gamma}{2})\right)^2}{\epsilon^2}$$
A: In the general case, the distribution generated by your deterministic function, based on the normally distributed inputs, may not itself be normal. Nor might it be solvable in closed form. However, in either or both cases, you can analyze the resulting distribution empirically, from the results of running the Monte Carlo distribution many times, and looking at the empirical variance, and other statistical measures of interest, of the resulting empirical samples.
Of course, if your distribution is analyzable in closed form, and is Gaussian, given Gaussian input, you might be able to derive a closed-form formulae for the variance of the output distribution, as a function of the variance and mean of the input Gaussian distribution.
A: There's no need to complicate this. Monitor the metric that you're calculating, then stop when it stops changing.
For instance, if you're calculating the probability of an event, then you monitor the quantity $$\hat p=\frac{n_i}{i},$$
where $n_i$ - number of observations of an event after $i$ simulations. Set a threshold such as 0.01 and then watch $\Delta p_i=p_i-p_{i-1}$ or $\delta p_i=\Delta p_i/p_i$, once $\Delta p_i<0.01$ or $\delta p_i<0.01$ stop the simulations. The former is the absolute error and the latter is relative error.
The final $i$ is the number you're looking for.
The only time this will not work if when you have to know $i$ in advance for some reason, e.g. for capacity planning purposes, and you absolutely cannot do any trial runs. In this case you have a problem. In the example above you could try approximations if you know the order of magnitude of $\hat p$. You could use Wald's formula for variance to back out $i$ like this:
$$\sigma=\sqrt{\hat p(1-\hat p)/i}$$
The problem's that if you're using variance reduction techniques (which you should) then these formulas don't work that well: they'll overestimate the number of simulations.
