function ZI = newinterp2 (X,Y,Z,XI,YI)
if (nargin == 3)
  [X, Y] = meshgrid(1:rows (Z), 1:columns (Z));
end
% X and Y are meshgrid(x,y)
% Z is a function of X and Y
X = X(1,:);
Y = Y(:,1);
nx = length (X);
ny = length (Y);
ppx = spline (X, eye (nx));
Ax = ppval (ppx, XI(:));
ppy = spline (Y, eye (ny));
Ay = ppval (ppy, YI(:));
if (size (XI, 1) == 1 && size (YI, 2) == 1)
% ZI corresponds to meshgrid (XI,YI)    
  ZI = Ay.' * Z * Ax;
else
% scattered 
  ZI = reshape (sum ((Ay' * Z) .* Ax', 2), size (XI));
end