In Line 253, when calculating the distances for more than two dimensions, the code was
d = sqrt(sum((X(i,:)-X(j,:)).^2));
I think it should be
d = sqrt(sum((X(i,:)-X(j,:)).^2,2));
as we need the sums of the rows but not the columns. After this correction, everything works perfectly in calculating 3-dimensional variograms.
Many thanks for your function. It is absolutely fabulous!