diff --git a/scripts/general/interp2.m b/scripts/general/interp2.m
--- a/scripts/general/interp2.m
+++ b/scripts/general/interp2.m
@@ -180,7 +180,12 @@
     error ("interp2: Y must be strictly monotonic");
   endif
 
-  if (any (strcmp (method, {"nearest", "linear", "pchip", "cubic"})))
+  if (strcmp (method, "cubic") && (rows (Z) < 3 || columns (Z) < 3))
+    warning ("interp2: cubic requires at least 3 points in each dimension. Falling back to linear interpolation.")
+    method = "linear";
+  endif
+
+  if (any (strcmp (method, {"nearest", "linear", "pchip"})))
 
     ## If Xi and Yi are vectors of different orientation build a grid
     if ((isrow (XI) && iscolumn (YI)) || (iscolumn (XI) && isrow (YI)))
@@ -244,11 +249,10 @@
       idx = sub2ind (size (Z), yidx+jj, xidx+ii);
       ZI = Z(idx);
 
-    elseif (strcmp (method, "pchip") || strcmp (method, "cubic"))
+    elseif (strcmp (method, "pchip"))
 
       if (length (X) < 2 || length (Y) < 2)
-        error ("interp2: %s requires at least 2 points in each dimension",
-               method);
+        error ("interp2: pchip requires at least 2 points in each dimension");
       endif
 
       ## first order derivatives
@@ -328,10 +332,45 @@
       else
         error ("interp2: XI, YI must have uniform spacing ('meshgrid' format)");
       endif
-    endif
-
-    return; # spline doesn't need NA extrapolation value (MATLAB compatibility)
-
+      return; # spline doesn't need NA extrapolation value (MATLAB compatibility)
+    elseif (strcmp (method, "cubic"))
+      ## check if interpolation points are a meshgrid and reduce to vectors if possible
+      if (size_equal (XI, YI) && all (all (XI(1, :) == XI & YI(:, 1) == YI)))
+        XI = XI(1, :);
+        YI = YI(:, 1);
+      endif
+
+      ## make X a row vector
+      X = X.';
+
+      ## quadratic padding + additional zeros for the special case of copying
+      ## the last line (like x=1:5, xi=5, requires to have indexes 6 and 7)
+      row_1 = 3*Z(1, :, :) - 3*Z(2, :, :) + Z(3, :, :);
+      row_end = 3*Z(end, :, :) - 3*Z(end-1, :, :) + Z(end-2, :, :);
+      ZI = [3*row_1(:, 1, :) - 3*row_1(:, 2, :) + row_1(:, 3, :), ...
+            row_1, ...
+            3*row_1(:, end, :) - 3*row_1(:, end-1, :) + row_1(:, end-2, :), ...
+            0;
+            #
+            3*Z(:, 1, :) - 3*Z(:, 2, :) + Z(:, 3, :), ...
+            Z, ...
+            3*Z(:, end, :) - 3*Z(:, end-1, :) + Z(:, end-2, :), ...
+            zeros(rows (Z), 1, size (Z, 3));
+            #
+            3*row_end(:, 1, :) - 3*row_end(:, 2, :) + row_end(:, 3, :), ...
+            row_end, ...
+            3*row_end(:, end, :) - 3*row_end(:, end-1, :) + row_end(:, end-2, :), ...
+            0;
+            zeros(1, columns (Z) + 3, size (Z, 3))];
+
+      ## interpolate
+      if (isrow (XI) && iscolumn (YI))
+        ZI = conv_interp_vec (ZI, Y, YI, @cubic, [-2, 2], 1);
+        ZI = conv_interp_vec (ZI, X, XI, @cubic, [-2, 2], 2);
+      else
+        ZI = conv_interp_pairs (ZI, X, Y, XI, YI, @cubic, [-2, 2]);
+      endif
+    endif
   endif
 
   ## extrapolation 'extrap'
@@ -364,6 +403,66 @@
   b = ! any (d1(:) != 0);
 endfunction
 
+## cubic convolution kernel with a = -0.5 for MATLAB compatibility.
+function w = cubic (h)
+  absh = abs (h);
+  absh01 = absh <= 1;
+  absh12 = absh <= 2 & ~absh01;
+  w = (1.5 * absh.^3 - 2.5 * absh.^2 + 1)             .* absh01 + ...  ## for |h| <= 1
+      (-0.5 * absh.^3 + 2.5 * absh.^2 - 4 * absh + 2) .* absh12;       ## for 1 < |h| <= 2
+end
+
+## bicubic interpolation of full matrix in one direction with vector
+function out = conv_interp_vec (Z, XY, XIYI, kernel, kernel_bounds, axis)
+  ## allocate output
+  out_shape = size (Z);
+  out_shape(axis) = length (XIYI);
+  out = zeros (out_shape);
+
+  ## get indexes and distances h
+  spread = abs (XY(1) - XY(2));
+  idx = lookup (XY, XIYI, "l");
+  h = (XIYI - XY(idx)) / spread;
+  idx += -kernel_bounds(1) - 1; ## apply padding for indexes
+
+  ## interpolate
+  for shift = kernel_bounds(1)+1:kernel_bounds(2)
+    if axis == 1
+      out += Z(idx + shift, :, :) .* kernel (shift - h);
+    else
+      out += Z(:, idx + shift, :) .* kernel (shift - h);
+    endif
+  endfor
+endfunction
+
+## bicubic interpolation of arbitrary XI-YI-pairs
+function out = conv_interp_pairs (Z, X, Y, XI, YI, kernel, kernel_bounds)
+  spread_x = abs (X(1, 1) - X(1, 2));
+  spread_y = abs (Y(1, 1) - Y(2, 1));
+  idx_x = lookup (X, XI, "l");
+  idx_y = lookup (Y, YI, "l");
+  h_x = (XI - reshape (X(idx_x), size (idx_x))) / spread_x;
+  h_y = (YI - reshape (Y(idx_y), size (idx_y))) / spread_y;
+
+  idx_x += -kernel_bounds(1) - 1; ## apply padding for indexes
+  idx_y += -kernel_bounds(1) - 1; ## apply padding for indexes
+
+  out = zeros (rows (YI), columns (XI), size (Z, 3));
+  for shift_x = kernel_bounds(1)+1:kernel_bounds(2)
+    ## interpolate in y-direction
+    out_y = zeros (rows (YI), columns (XI), size (Z, 3));
+    for shift_y = kernel_bounds(1)+1:kernel_bounds(2)
+      pixels = Z(sub2ind (size (Z), idx_y + shift_y, idx_x + shift_x));
+      out_y += pixels .* kernel(shift_y - h_y);
+    endfor
+
+    ## interpolate in x-direction
+    out += out_y .* kernel (shift_x - h_x);
+  endfor
+endfunction
+
+
+
 
 %!demo
 %! clf;
@@ -409,8 +508,7 @@
 %! [x,y] = meshgrid (x,y);
 %! hold on; plot3 (x,y,A,"b*"); hold off;
 
-## 'pchip' commented out since it is the same as 'cubic'
-%!#demo
+%!demo
 %! clf;
 %! colormap ("default");
 %! A = [13,-1,12;5,4,3;1,6,2];
@@ -421,8 +519,7 @@
 %! [x,y] = meshgrid (x,y);
 %! hold on; plot3 (x,y,A,"b*"); hold off;
 
-## 'pchip' commented out since it is the same as 'cubic'
-%!#demo
+%!demo
 %! clf;
 %! colormap ("default");
 %! [x,y,A] = peaks (10);
@@ -477,21 +574,100 @@
 %! [x,y] = meshgrid (x,y);
 %! hold on; plot3 (x,y,A,"b*"); hold off;
 
+%!shared x, y, orig, xi, yi, expected
 %!test  # simple test
 %! x = [1,2,3];
 %! y = [4,5,6,7];
 %! [X, Y] = meshgrid (x, y);
 %! orig = X.^2 + Y.^3;
-%! xi = [1.2,2, 1.5];
-%! yi = [6.2, 4.0, 5.0]';
-%!
+%! xi = [0.8, 1.2, 2.0, 1.5];
+%! yi = [6.2, 4.0, 5.0, 7.1]';
+%!
+%! # check nearest neighbor
 %! expected = ...
-%!   [243,   245.4,  243.9;
-%!     65.6,  68,     66.5;
-%!    126.6, 129,    127.5];
-%! result = interp2 (x,y,orig, xi, yi);
-%!
+%!   [NA, 217, 220, 220;
+%!    NA,  65,  68,  68;
+%!    NA, 126, 129, 129;
+%!    NA,  NA,  NA,  NA];
+%! result = interp2 (x, y, orig, xi, yi, "nearest");
+%! assert (result, expected);
+%!
+%! # check invariance of translation
+%! result = interp2 (x+3, y-7, orig, xi+3, yi-7, "nearest");
+%! assert (result, expected);
+%!
+%! # check invariance of scaling
+%! result = interp2 (x*3, y*(-7), orig, xi*3, yi*(-7), "nearest");
+%! assert (result, expected);
+%!
+%! # check interpolation with index pairs
+%! result = interp2 (x, y, orig, xi(2:4), yi(1:3)', "nearest");
+%! assert (result, expected(sub2ind(size(expected), 1:3, 2:4)));
+%!
+%! # check bilinear interpolation
+%! expected = ...
+%!   [NA, 243,     245.4,   243.9;
+%!    NA,  65.6,    68,      66.5;
+%!    NA, 126.6,   129,     127.5;
+%!    NA,  NA,      NA,      NA];
+%! result = interp2 (x, y, orig, xi, yi);
 %! assert (result, expected, 1000*eps);
+%!
+%! # check invariance of translation
+%! result = interp2 (x+3, y-7, orig, xi+3, yi-7);
+%! assert (result, expected, 1000*eps);
+%!
+%! # check invariance of scaling
+%! result = interp2 (x*3, y*(-7), orig, xi*3, yi*(-7));
+%! assert (result, expected, 1000*eps);
+%!
+%! # check interpolation with index pairs
+%! result = interp2 (x, y, orig, xi(2:4), yi(1:3)');
+%! assert (result, expected(sub2ind(size(expected), 1:3, 2:4)), 1000*eps);
+%!
+%! # check spline interpolation
+%! expected = ...
+%!   [238.968,  239.768,  242.328,  240.578;
+%!     64.64,    65.44,    68,       66.25;
+%!    125.64,   126.44,   129,      127.25;
+%!    358.551,  359.351,  361.911,  360.161];
+%! result = interp2 (x, y, orig, xi, yi, "spline");
+%! assert (result, expected, 1000*eps);
+%!
+%! # check invariance of translation
+%! result = interp2 (x+3, y-7, orig, xi+3, yi-7, "spline");
+%! assert (result, expected, 1000*eps);
+%!
+%! # check invariance of scaling
+%! result = interp2 (x*3, y*(-7), orig, xi*3, yi*(-7), "spline");
+%! assert (result, expected, 1000*eps);
+%!
+%!xtest
+%! # FIXME: spline interpolation does not support index pairs, Matlab does.
+%! result = interp2 (x, y, orig, xi(2:4), yi(1:3)', "spline");
+%! assert (result, expected(sub2ind(size(expected), 1:3, 2:4)), 1000*eps);
+%!
+%!test
+%! # check bicubic interpolation
+%! expected = ...
+%!   [NA, 239.96,  242.52,  240.77;
+%!    NA,  65.44,   68,      66.25;
+%!    NA, 126.44,  129,     127.25;
+%!    NA,  NA,      NA,      NA];
+%! result = interp2 (x, y, orig, xi, yi, "cubic");
+%! assert (result, expected, 10000*eps);
+%!
+%! # check invariance of translation
+%! result = interp2 (x+3, y-7, orig, xi+3, yi-7, "cubic");
+%! assert (result, expected, 10000*eps);
+%!
+%! # check invariance of scaling
+%! result = interp2 (x*3, y*(-7), orig, xi*3, yi*(-7), "cubic");
+%! assert (result, expected, 10000*eps);
+%!
+%! # check interpolation with index pairs
+%! result = interp2 (x, y, orig, xi(2:4), yi(1:3)', "cubic");
+%! assert (result, expected(sub2ind(size(expected), 1:3, 2:4)), 10000*eps);
 
 ## Test that interpolating a complex matrix is equivalent to interpolating its
 ## real and imaginary parts separately.
@@ -499,8 +675,8 @@
 %!
 %! xi = [2.5, 3.5];
 %! yi = [0.5, 1.5]';
-%! orig = rand (2, 3) + i * rand (2, 3);
-%! for method = {"nearest", "linear", "pchip", "spline"}
+%! orig = rand (4, 3) + i * rand (4, 3);
+%! for method = {"nearest", "linear", "pchip", "cubic", "spline"}
 %!   interp_complex = interp2 (orig, xi, yi, method{1});
 %!   interp_real = interp2 (real (orig), xi, yi, method{1});
 %!   interp_imag = interp2 (imag (orig), xi, yi, method{1});
@@ -540,6 +716,8 @@
 %! assert (interp2 (x,y,orig, xi, yi,"spline", 2), [2,2;2,2]);
 %! assert (interp2 (x,y,orig, xi, yi,"linear", 0+1i), [0+1i,0+1i;0+1i,0+1i]);
 %! assert (interp2 (x,y,orig, xi, yi,"spline"), [27,43;512,528]);
+%! assert (interp2 (x,y,orig, xi, yi,"cubic"), [NA,NA;NA,NA]);
+%! assert (interp2 (x,y,orig, xi, yi,"cubic", 2), [2,2;2,2]);
 
 %!test  # for values at boundaries
 %! A = [1,2;3,4];
@@ -548,7 +726,7 @@
 %! assert (interp2 (x,y,A,x,y,"linear"), A);
 %! assert (interp2 (x,y,A,x,y,"nearest"), A);
 
-%!test  # for Matlab-compatible rounding for 'nearest'
+%!test  # for Matlab-compatible rounding for "nearest"
 %! X = meshgrid (1:4);
 %! assert (interp2 (X, 2.5, 2.5, "nearest"), 3);
 
@@ -567,20 +745,26 @@
 %!assert (interp2 (z), zout, tol)
 %!assert (interp2 (z, "linear"), zout, tol)
 %!assert (interp2 (z, "pchip"), zout, tol)
-%!assert (interp2 (z, "cubic"), zout, 10 * tol)
+%!assert (interp2 (z, "cubic"), zout, tol)
 %!assert (interp2 (z, "spline"), zout, tol)
 %!assert (interp2 (z, [2 3 1], [2 2 2]', "linear"),
 %!        repmat ([5, 7, 3], [3, 1]), tol)
 %!assert (interp2 (z, [2 3 1], [2 2 2]', "pchip"),
 %!        repmat ([5, 7, 3], [3, 1]), tol)
 %!assert (interp2 (z, [2 3 1], [2 2 2]', "cubic"),
-%!        repmat ([5, 7, 3], [3, 1]), 10 * tol)
+%!        repmat ([5, 7, 3], [3, 1]), tol)
 %!assert (interp2 (z, [2 3 1], [2 2 2]', "spline"),
 %!        repmat ([5, 7, 3], [3, 1]), tol)
 %!assert (interp2 (z, [2 3 1], [2 2 2], "linear"), [5 7 3], tol)
 %!assert (interp2 (z, [2 3 1], [2 2 2], "pchip"), [5 7 3], tol)
-%!assert (interp2 (z, [2 3 1], [2 2 2], "cubic"), [5 7 3], 10 * tol)
+%!assert (interp2 (z, [2 3 1], [2 2 2], "cubic"), [5 7 3], tol)
 %!assert (interp2 (z, [2 3 1], [2 2 2], "spline"), [5 7 3], tol)
+%!assert (interp2 (z, [3; 3; 3], [2; 3; 1], "linear"), [7; 9; 5], tol)
+%!assert (interp2 (z, [3; 3; 3], [2; 3; 1], "pchip"), [7; 9; 5], tol)
+%!assert (interp2 (z, [3; 3; 3], [2; 3; 1], "cubic"), [7; 9; 5], tol)
+%!xtest
+%! # FIXME: single column yields single row with spline interpolation (numbers are correct)
+%! assert (interp2 (z, [3; 3; 3], [2; 3; 1], "spline"), [7; 9; 5], tol)
 
 ## Test input validation
 %!error interp2 (1, 1, 1, 1, 1, 2)    # only 5 numeric inputs
@@ -592,10 +776,10 @@
 %!error <N must be an integer .= 0> interp2 (ones (2), -1)
 %!error <N must be an integer .= 0> interp2 (ones (2), 1.5)
 %!warning <ignoring unsupported '\*' flag> interp2 (rand (3,3), 1, "*linear");
-%!error <EXTRAP must be a numeric scalar> interp2 (1, 1, 1, 1, 1, 'linear', {1})
-%!error <EXTRAP must be a numeric scalar> interp2 (1, 1, 1, 1, 1, 'linear', ones (2,2))
-%!error <EXTRAP must be a numeric scalar> interp2 (1, 1, 1, 1, 1, 'linear', "abc")
-%!error <EXTRAP must be a numeric scalar> interp2 (1, 1, 1, 1, 1, 'linear', "extrap")
+%!error <EXTRAP must be a numeric scalar> interp2 (1, 1, 1, 1, 1, "linear", {1})
+%!error <EXTRAP must be a numeric scalar> interp2 (1, 1, 1, 1, 1, "linear", ones (2,2))
+%!error <EXTRAP must be a numeric scalar> interp2 (1, 1, 1, 1, 1, "linear", "abc")
+%!error <EXTRAP must be a numeric scalar> interp2 (1, 1, 1, 1, 1, "linear", "extrap")
 %!error <X, Y must be numeric matrices> interp2 ({1}, 1, ones (2), 1, 1)
 %!error <X, Y must be numeric matrices> interp2 (1, {1}, ones (2), 1, 1)
 %!error <XI, YI must be numeric> interp2 (1, 1, ones (2), {1}, 1)
@@ -611,3 +795,4 @@
 %!error <XI, YI must have uniform spacing> interp2 (1:2, 1:2, ones (2), [1 2 4], [1 2 3], "spline")
 %!error <XI, YI must have uniform spacing> interp2 (1:2, 1:2, ones (2), [1 2 3], [1 2 4], "spline")
 %!error interp2 (1, 1, 1, 1, 1, "foobar")
+%!warning <cubic requires at least 3 points in each dimension.> interp2 (eye(2), 1.5, 1.5, "cubic");

