@Wolfgang: I really appreciate your interactivity and your pertinence in all your comments. I just coded the Deutsch method to solve the problem of the negative weights. Thanks to your doi link. Thus, the gaussian variogram is now well working :)
Here are the lines code to insert within the "for" loop, between the paragraphs "solve system" and "estimate zi".
However, I'm not an excellent developper. I'll let you find some better/faster way to code this.
% solve system
lambda = A*b;
% correct negative weights (Deutsch 1996)
if sum(lambda(:)<0) % if negative weigths are found
% compute the "average absolute magnitude of the negative weights"
lambdaLim=lambda(1:end-1,:);
PosInd=lambdaLim>=0;
lambdaLim(PosInd)=NaN;
lambdaLim=nanmean(abs(lambdaLim),1);
% compute the mean covariance between "negative" locations and the IX ones
bNoOnes=b(1:end-1,:);
bLim=bNoOnes;
bLim(PosInd)=NaN;
bLim=nanmean(bLim,1);
% correct the weights
lambdaC=lambda(1:end-1,:); % step 1 of Deutsch 1996
lambdaC(lambdaC<0)=0; % step 2 of Deutsch 1996
lambdaC(PosInd & bNoOnes-repmat(bLim,numobs,1)<0 & lambda(1:end-1,:)-repmat(lambdaLim,numobs,1)<0)=0; % step 3 of Deutsch 1996
% restandardization of the weights
lambda(1:end-1,:)=lambdaC./repmat(sum(lambdaC,1),numobs,1);
end
% estimate zi
zi(IX) = lambda'*z;

@Sébastien: Thanks for pointing at the error when using a nugget variogram.
I am currently abroad, so I probably won't find much time to update the function. So, in the meantime please adopt the changes proposed by Sébastien.
The problem with the gaussian variogram is most likely due to negative weights (see http://dx.doi.org/10.1016/0098-3004(96)00005-2 ). I will try to implement the procedure suggested by Deutsch (1996) in a future update.

Very nice submission. I was amazed by its speed.
However, two things to say:
1) there is a mistake in the code when expanding b with ones. It should be the following:
% expand b with ones
b = vstruct.func([vstruct.range vstruct.sill],b);
if ~isempty(vstruct.nugget)
b = b+vstruct.nugget;
end
b=[b;ones(1,chunksize)];
Indeed, the ones should be added AFTER adding the nugget. Otherwise, the nugget is also applied to the ones: an offset is then unfortunately observed on the results...
2) I have the same problem than Kirsten: it's working well with spherical variograms. However, when I apply a gaussian one, huge values are found. I'm actually investigating the problem. If I find something, I'll let you know. But I would appreciate if you find the solution first ;)

@Imran Zafar: Fig. 4 shows the kriging variance which quantifies the uncertainty of the kriging estimates. You will note, that the kriging variance increases with increasing distance from observations.

Comment only