Let $X_i := B(t_i) [B(t_{i+1}) - B(t_i)]$ for $i=0$ to $n-1$. The rv.s
$X_i$ are gaussian and independent due to the independence of the
increments of $B(t)$. Moreover $\mathbb{E}(X_i) = 0$ and
$$
\text{Var}(X_i) =
\mathbb{E}[B(t_i)] \, \mathbb{E}[B(t_{i+1}) - B(t_i)] = t_i \, (t_{i+1} - t_i),
%% = \frac{i}{n} \frac{1}{n} = \frac{i}{n^2}
$$$$
\text{Var}(X_i) = E[X_i^2] =
\mathbb{E}[B(t_i)^2] \, \mathbb{E}[(B(t_{i+1}) - B(t_i))^2] = t_i \, (t_{i+1} - t_i),
%% = \frac{i}{n} \frac{1}{n} = \frac{i}{n^2}
$$
so $\text{Var}(X_i) = i / n^2$ if the points $t_i := i /n$ are used. Then
$$
\text{Var}(I^{(n)}) = \sum_{i=0}^{n-1} \text{Var}(X_i) =
\frac{1}{n^2} \sum_{i=0}^{n-1} i
= \frac{1}{n^2} \, \frac{(n-1)n}{2} = \left[1 - \frac{1}{n}\right] \frac{1}{2}.
$$
Up to the term $1/n$, the distribution of $I^{(n)}$ remains nearly the
same for all $n$. By drawing a number, say $N_{\text{sim}}$, of
independent $I^{(n)}$, we get a sample of size $N_{\text{sim}}$ from a
distribution depending only slightly of $n$ and which is nearly the
normal with mean $0$ and variance $1/2$. So in the program, the
elements of v_sim corresponding to n are (unless n is small)
simply the variances of independent samples of size $N_{\text{sim}}$
from $\text{Norm}(0, \,1/2)$.
The number $n$ of points is nearly irrelevant here because of the specific integrand considered. The situation can be compared to the evaluation of the non-stochastic integral $\int_0^1 t \,\text{d}t$ using the trapezoidal rule. The result will not change when more trapezes are used. Obviously things would be different if a different adapted process or function was used as integrand, a denser design $[t_i]_i$ leading then to a smaller bias.
