ordinary kriging power method

i'm trying to generate simple coding for ordinary kriging method. this method is using using power variogram model (power=1.5). This is use for estimate unknown value of z at (x,y) by using known (x,y,z values). however the results seem to be not accurate as the interpolation value given at known point is not the same as the original one. For example, at (1,2) the value is 7. The return calculation give the answer at (1,2) z=4.369. Maybe there is something missing/wrong in the looping process? This is the code i'm working on right now. Thank you.

%import data from .txt file
Pressure=importdata('xyz.txt');
Coordinate=Pressure.data;
x=Coordinate(:,1);
y=Coordinate(:,2);
z=Coordinate(:,3);

%range of unknown value coordinates
%predicitons matrix
[X,Y] = meshgrid(linspace(0,10,11),linspace(0,10,11));
xi=X(:);
yi=Y(:);

% size of input arguments
xi = xi(:);
yi = yi(:);
numest = numel(xi);
numobs = numel(x);
chunksize=numest;

% distance matrix of locations with known values
Dx = hypot(bsxfun(@minus,x,x'),bsxfun(@minus,y,y'));
%Data covariance matrix using power variogram model (power=1.5)
A=Dx.^(1.5);

% matrix expanded by one line and one row to for ordinary kriging
% condition, that all weights must sum to one (lagrange multiplier)
A = [[A ones(numobs,1)];ones(1,numobs) 0];

% expand z
z  = [z;0];

%output zi
zi = nan(numest,1);
nrloops   = ceil(numest/chunksize);

% looping 
for r = 1:nrloops;
    % built chunks
        IX = (r-1)*chunksize +1 : numest;
        chunksize = numel(IX);
end
 b = hypot(bsxfun(@minus,x,xi(IX)'),bsxfun(@minus,y,yi(IX)'));
 
 % expand b with ones
b = [b;ones(1,chunksize)];

%calculate lambda
lambda = A^-1*b;

%calculate z (predicted) values
zi(IX)  = lambda'*z;

%sort into xyz coordinate
Z=[xi yi zi];