Using a for loop to manually construct the correlation matrix and a bivariate normal distribution is likely very slow. In general, to speed up the computation that involves multivariate normal, it is necessary to find a way to express the covariant matrix as using its Cholesky decomposition.
There are a few similar tricks in GP, you might find posts here helpful such as: Multiple (uncertain) function observations of the same Gaussian process