Misconception about Long-Run Variance AR1 Matlab

Question

I think I have some misconception about long-run variances in AR(1) processes, maybe you can help out. I am simulating an AR1 process with persistance 0.5 using the ARMA function in Matlab over a thousand sample paths and 1000 obs each.

ar = arima('Constant',c,'AR',{0.5},'Variance',1);
Y = simulate(ar,1000, 'NumPaths',1000);

In my Advanced Econometrics class, I learned that one can calculate the long-run variable of a covariance-stationary process by summing up the auto-covariances between the -infinite and +infinite order. Calculating the covariances I rely on xcorr such that the below snippet generates all auto-covariances for all available lags for each reps individually (thus creates a T*2 x R Matrix). According to the idea of summing auto-covariances, I calculate a cumulative sum of covariances too.

for r = 1:R
    covY(:,r) = xcorr(Y(:,r), 'biased');
    sum_covY = cumsum(covY,1);
end

I plot this cumulative sum of covariancs for increasing lags for all R sample paths and also plot the average across all sample paths accordingly.

plot(-2*N/2+1:2*N/2-1, sum_covY,'Color',[.85,.85,.85]), xlabel('Lags', 'Interpreter','latex'), title('Long-Run Variance','Interpreter','latex')
hold on
yline(sigma^2/(1-(rho^2)))
h = plot(-2*N/2+1:2*N/2-1, mean(sum_covY,2),'k','LineWidth',2);
legend(h,'Simulation Mean','Location','NorthWest');
hold off 

Now, the misconception that I have. Should this Long-Run Variance not converge to what calculating the variance based on the AR(1) equation, say, y_t = 0.5*y_t-1 + eps should yield? (With eps ~ N(0,1)). That would be Cov(y_t, y_t-1) = sigma^2 / 1 - 0.5^2 = 1.333. Can someone please shed some light on this? The black line in the attached image converges to 4.15.

Thank you so much!

Marcello