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);
b = b+vstruct.nugget;
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 ;)