@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