function ZI = newinterp2 (varargin)
if (nargin == 3)
  Z = varargin{1};
  XI = varargin{2};
  YI = varargin{3};
  X = 1:rows(Z);
  Y = 1:columns (Z);
  method = "cubic";
elseif (nargin == 4)
  Z = varargin{1};
  XI = varargin{2};
  YI = varargin{3};
  X = 1:rows(Z);
  Y = 1:columns (Z);
  method = varargin{4};
elseif (nargin == 5)
  X = varargin{1};
  Y = varargin{2};
  Z = varargin{3};
## X is a row vector, Y is a column vector 
## Z corresponds to meshgrid (X,Y)
  XI = varargin{4};
  YI = varargin{5};
  method = "cubic";
elseif (nargin == 6)
  X = varargin{1};
  Y = varargin{2};
  Z = varargin{3};
## X is a row vector, Y is a column vector 
## Z corresponds to meshgrid (X,Y)
  XI = varargin{4};
  YI = varargin{5};
  method = varargin{6};
endif
Y = Y.';
if (any (abs (diff (X) - diff (X)(1)) > eps) || ...
    any (abs (diff (Y) - diff (Y)(1)) > eps) || ...
    strcmp (method, "spline"))
## not equispaced X, Y or spline required: spline interpolation
  if (all (abs (diff (XI) - diff (XI)(1)) < eps) && ...
    all (abs (diff (YI) - diff (YI)(1)) < eps))
## equispaced XI, YI: __splinen__
    ZI = __splinen__ ({Y, X},Z,{YI, XI},0,"spline");
  else
## not equispaced XI, YI: spline
    Ax = spline (X, eye (length (X)), XI(:));
    Ay = spline (Y, eye (length (Y)), YI(:));
    ZI = reshape (sum ((Ay' * Z) .* Ax', 2), size (XI));
  endif
else
## equispaced X, Y and spline not required: bicubic interpolation
  ZI = newbicubic (X, Y, Z, XI, YI);
endif
endfunction