diff --git a/inst/imremap.m b/inst/imremap.m
--- a/inst/imremap.m
+++ b/inst/imremap.m
@@ -16,6 +16,7 @@
 ## -*- texinfo -*-
 ## @deftypefn {Function File} @var{warped} = imremap(@var{im}, @var{XI}, @var{YI})
 ## @deftypefnx{Function File} @var{warped} = imremap(@var{im}, @var{XI}, @var{YI}, @var{interp}, @var{extrapval})
+## @deftypefnx{Function File} @var{warped} = imremap(@var{im}, @var{XI}, @var{YI}, 'bicubic', @var{padding})
 ## Applies any geometric transformation to the image @var{im}.
 ##
 ## The arguments @var{XI} and @var{YI} are lookup tables that define the resulting
@@ -37,64 +38,191 @@
 ##
 ## All values of the result that fall outside the original image will
 ## be set to @var{extrapval}.  The default value of @var{extrapval} is 0.
+## For bicubic interpolation it is possible to apply @var{padding} instead.
+## Valid padding methods are: 'replicate', 'symmetric', 'reflect', 'circular'.
 ##
 ## @seealso{imperspectivewarp, imrotate, imresize, imshear, interp2}
 ## @end deftypefn
 
 function [warped] = imremap(im, XI, YI, interp = "linear", extrapval = 0)
+  interp = interp_method (interp);
 
   if (nargin < 3 || nargin > 5)
     print_usage ();
   elseif (!  ((isnumeric (im) || islogical (im)) && ! issparse (im)
             && ! isempty (im)) || ndims (im) > 3)
     error ("imremap: IM must be a grayscale or RGB image.")
-  elseif (! size_equal (XI, YI) || ! ismatrix (XI) || ! isnumeric (XI))
-    error ("imremap: XI and YI must be matrices of the same size");
+  elseif (! (size_equal (XI, YI) || (isvector (XI) && isvector (YI))) || ! ismatrix (XI) || ! isnumeric (XI))
+    error ("imremap: XI and YI must be matrices of the same size or vectors");
   elseif (! ischar (interp))
     error ("imremap: INTERP must be a string with interpolation method")
-  elseif (! isscalar (extrapval))
-    error ("imremap: EXTRAPVAL must be a scalar");
+  elseif (! isscalar (extrapval) && ! (ischar (extrapval) && strcmp (interp, 'cubic')))
+    error ("imremap: Specify a scalar EXTRAPVAL for constant padding or in case of bicubic interpolation a string for the PADDING method");
+  endif
+
+  ## check if interpolation points are a meshgrid and reduce to vectors if possible
+  if (! isvector (XI) && size_equal (XI, YI) && all (all (repmat (XI(1, :), [rows(XI), 1]) == XI & repmat (YI(:, 1), [1, columns(YI)]) == YI)))
+    XI = XI(1, :);
+    YI = YI(:, 1);
+  endif
+
+  ## if XI and YI are vectors, make sure XI is a row vector and YI a column vector for broadcasting
+  if (iscolumn (XI) && !isvector (im))
+    XI = XI';
   endif
-  interp = interp_method (interp);
+  if (isrow (YI) && !isvector (im))
+    YI = YI';
+  endif
+
+  ## for bicubic interpolation do not use interp2, but another implementation for Matlab compatibility
+  if (strcmp (interp, "cubic"))
+    padding = 'symmetric';
+    if (ischar (extrapval))
+      padding = extrapval;
+    endif
+
+    ## interpolate
+    warped = bicubic_conv (double (im), XI, YI, padding);
 
-  ## Interpolate
-  sz = size (im);
-  n_planes = prod (sz(3:end));
-  sz(1:2) = size (XI);
-  warped = zeros(sz);
-  for i = 1:n_planes
-    # 1-pixel image planes:
-    if isscalar (im(:,:,i))
-      if (all (XI(:)==1) && all (YI(:)==1))
-        warped(:,:,i) = double (im(:,:,i));
+    if (isscalar (extrapval))
+      ## values got padded for smooth borders, but constant padding has been requested
+      outside = (XI < 0.5) | (XI > columns (im) + 0.5) | ...
+                (YI < 0.5) | (YI > rows (im) + 0.5);
+      outside = repmat (outside, [1, 1, size(im,3)]);
+      warped(outside) = extrapval;
+    endif
+
+  else
+    sz = size (im);
+    n_planes = prod (sz(3:end));
+    sz(1) = size (YI, 1);
+    sz(2) = size (XI, 2);
+    warped = zeros (sz);
+    for i = 1:n_planes
+      # 1-pixel image planes:
+      if isscalar (im(:,:,i))
+        if (all (XI(:)==1) && all (YI(:)==1))
+          warped(:,:,i) = double (im(:,:,i));
+        else
+          warped(:,:,i) = extrapval;
+        endif
+      # row image planes:
+      elseif size(im, 1) == 1
+        if (all (YI(:)==1))
+          warped(:,:,i) = ones (sz(1), 1) * interp1 ([1:size(im,2)], double (im(1,:,i)), XI, interp, extrapval);
+        else
+          warped(:,:,i) = extrapval;
+        endif
+      # col image planes:
+      elseif size(im, 2) == 1
+        if (all (XI(:)==1))
+          warped(:,:,i) = interp1 ([1:size(im,1)], double (im(:,1,i)), YI, interp, extrapval) * ones (1, sz(2));
+        else
+          warped(:,:,i) = extrapval;
+        endif
+      # 2d image planes:
       else
-        warped(:,:,i) = extrapval;
+        warped(:,:,i) = interp2 (double (im(:,:,i)), XI, YI, interp, extrapval);
       endif
-    # line image planes:
-    elseif size(im, 1) == 1
-      if (all (YI(:)==1))
-        warped(:,:,i) = ones (size(XI,1), 1) * interp1 ([1:size(im,2)], double (im(1,:,i)), XI(1,:), interp, extrapval);
-      else
-        warped(:,:,i) = extrapval;
-      endif
-    # row image planes:
-    elseif size(im, 2) == 1
-      if (all (XI(:)==1))
-        warped(:,:,i) = interp1 ([1:size(im,1)], double (im(:,1,i)), YI(:,1), interp, extrapval) * ones (1, size(YI,2));
-      else
-        warped(:,:,i) = extrapval;
-      endif
-    # 2d image planes:
-    else
-      warped(:,:,i) = interp2 (double (im(:,:,i)), XI, YI, interp, extrapval);
-    endif
-  endfor
+    endfor
+  endif
 
   ## we return image on same class as input
   warped = cast (warped, class (im));
+endfunction
 
+## Cubic interpolation in 1d using a convolution kernel with a = -0.5 for MATLAB compatibility. Bicubic interpolation in interp2 is not MATLAB compatible.
+function w = cubic01(d, a)
+  ## requires: all (abs (d(:)) <= 1)
+  absd = abs (d);
+  w = (a+2) * absd.^3 - (a+3) * absd.^2 + 1;
+endfunction
+
+function w = cubic12(d, a)
+  ## requires all (1 < abs (d(:)) && abs (d(:)) <= 2)
+  absd = abs (d);
+  w = a * absd.^3 - 5*a * absd.^2 + 8*a * absd - 4*a;
+endfunction
+
+function p = intpolcub(I1, I2, I3, I4, D, a = -0.5)
+  ## requires: all (0 <= D(:) && D(:) < 1) && a < 1 && size_equal(I1, I2) && size_equal(I1, I3) && size_equal(I1, I4)
+  p = I1 .* cubic12 (-1-D, a) + I2 .* cubic01 (-D, a) + I3 .* cubic01 (1-D, a) + I4 .* cubic12 (2-D, a);
+endfunction
+
+## padding by changing indices. Cannot mimic constant value padding, like zero padding
+function idx = pad_indices(i, sz, method = 'symmetric')
+  if strcmp (method, 'replicate')
+    idx = max (min (i, sz), 1);
+  elseif strcmp (method, 'symmetric')
+    idx = i - 1;
+    m = mod (idx, sz);
+    odd = mod (floor (idx / sz), 2) == 1;
+    idx(odd) = sz - m(odd);
+    idx(!odd) = m(!odd) + 1;
+  elseif strcmp (method, 'reflect')
+    idx = i - 1;
+    while (any (idx(:) < 0 | idx(:) >= sz))
+      idx(idx < 0) = -idx(idx < 0);
+      idx(idx >= sz) = 2*sz - 2 - idx(idx >= sz);
+    endwhile
+    idx += 1;
+  elseif strcmp (method, 'circular')
+    idx = mod (i - 1, sz) + 1;
+  else
+    error (['Invalid argument for PADDING. Valid are "replicate", "symmetric", "reflect", "circular". You gave "', method, '"'])
+  endif
 endfunction
 
+## extract elements using meshgrid from 2d or 3d matrix
+function B = idx3(A, Y, X)
+  sz = size(A);
+  if (length (sz) == 2)
+    i = sub2ind (sz, Y, X);
+  else
+    X3 = repmat (X, [1, 1, sz(3)]);
+    Y3 = repmat (Y, [1, 1, sz(3)]);
+    ## Z is simply a 3d matrix, where the first channel is full of 1s, the second full of 2s, etc.
+    Z3 = reshape (repelem (1:sz(3), size (X, 1) * size (X, 2), 1), [size(X,1), size(X,2), sz(3)]);
+    i = sub2ind (sz, Y3, X3, Z3);
+  endif
+  B = A(i);
+endfunction
+
+## bicubic interpolation using convolution kernel
+function out = bicubic_conv(img, XI, YI, padding = 'symmetric')
+  ## make padded indices for smooth borders
+  K = floor(XI);
+  DX = XI - K;
+  K_1 = pad_indices (K - 1, columns (img), padding);
+  K0  = pad_indices (K    , columns (img), padding);
+  K1  = pad_indices (K + 1, columns (img), padding);
+  K2  = pad_indices (K + 2, columns (img), padding);
+
+  L = floor(YI);
+  DY = YI - L;
+  L_1 = pad_indices (L - 1, rows (img), padding);
+  L0  = pad_indices (L    , rows (img), padding);
+  L1  = pad_indices (L + 1, rows (img), padding);
+  L2  = pad_indices (L + 2, rows (img), padding);
+
+  if (isvector (XI) && !isvector (img)) ## rectilinear interpolation grid using vectors
+    ## interpolate in y-direction
+    out_y = intpolcub (img(L_1, :, :), img(L0, :, :), img(L1, :, :), img(L2, :, :), DY);
+
+    ## interpolate in x-direction
+    out = intpolcub (out_y(:, K_1, :), out_y(:, K0, :), out_y(:, K1, :), out_y(:, K2, :), DX);
+
+  else ## meshgrid interpolation
+    ## interpolate in y-direction at gridpoints
+    out_y_1 = intpolcub (idx3 (img, L_1, K_1), idx3 (img, L0, K_1), idx3 (img, L1, K_1), idx3 (img, L2, K_1), DY);
+    out_y0  = intpolcub (idx3 (img, L_1, K0),  idx3 (img, L0, K0),  idx3 (img, L1, K0),  idx3 (img, L2, K0),  DY);
+    out_y1  = intpolcub (idx3 (img, L_1, K1),  idx3 (img, L0, K1),  idx3 (img, L1, K1),  idx3 (img, L2, K1),  DY);
+    out_y2  = intpolcub (idx3 (img, L_1, K2),  idx3 (img, L0, K2),  idx3 (img, L1, K2),  idx3 (img, L2, K2),  DY);
+
+    ## interpolate in x-direction
+    out = intpolcub (out_y_1, out_y0, out_y1, out_y2, DX);
+  endif
+endfunction
 
 %!demo
 %! ## Generate a synthetic image and show it
@@ -120,3 +248,28 @@ endfunction
 %! XI = reshape(RXY(:,1), [100, 100]); YI = reshape(RXY(:,2), [100, 100]);
 %! warped = imremap(I, XI, YI);
 %! figure, imshow(warped);
+
+%!test
+%!
+## Test padding indirectly
+%! I = repmat([                      1  2  3  4                     ], [4, 1]);
+%! xi       = [-6 -5 -4 -3 -2 -1  0  1  2  3  4  5  6  7  8  9 10 11];
+%! exp_rep  = [ 1  1  1  1  1  1  1  1  2  3  4  4  4  4  4  4  4  4];
+%! exp_sym  = [ 2  3  4  4  3  2  1  1  2  3  4  4  3  2  1  1  2  3];
+%! exp_ref  = [ 2  1  2  3  4  3  2  1  2  3  4  3  2  1  2  3  4  3];
+%! exp_cir  = [ 2  3  4  1  2  3  4  1  2  3  4  1  2  3  4  1  2  3];
+%! yi       = 2.5;
+%!
+%! # rectilinear grid codepath
+%! assert (imremap (I, xi, yi, 'bicubic', 'replicate'), exp_rep);
+%! assert (imremap (I, xi, yi, 'bicubic', 'symmetric'), exp_sym);
+%! assert (imremap (I, xi, yi, 'bicubic', 'reflect'),   exp_ref);
+%! assert (imremap (I, xi, yi, 'bicubic', 'circular'),  exp_cir);
+%!
+%! # meshgrid codepath
+%! XI = [xi/2; xi; xi/2]; % cannot be reduced to vector, we will assert only middle row
+%! YI = repmat ([1.5; yi; 3.5], [1, length(xi)]);
+%! assert (imremap (I, XI, YI, 'bicubic', 'replicate')(2,:), exp_rep);
+%! assert (imremap (I, XI, YI, 'bicubic', 'symmetric')(2,:), exp_sym);
+%! assert (imremap (I, XI, YI, 'bicubic', 'reflect')(2,:),   exp_ref);
+%! assert (imremap (I, XI, YI, 'bicubic', 'circular')(2,:),  exp_cir);
diff --git a/inst/imresize.m b/inst/imresize.m
--- a/inst/imresize.m
+++ b/inst/imresize.m
@@ -104,20 +104,24 @@ function im = imresize (im, scale, metho
   idx_rows = 0.5 + off_rows + (0:out_rows-1) / scale_rows;
   idx_cols = 0.5 + off_cols + (0:out_cols-1) / scale_cols;
 
-  ## work around for imremap, valid for linear interpolation, probably bad
-  ## for cubic interpolation. Limits indices to the bounds 1 and end.
-  idx_rows = max (min (idx_rows, in_rows), 1);
-  idx_cols = max (min (idx_cols, in_cols), 1);
+  ## Simple replicate padding for nearest neighbor and bilinear interpolation. Padding for bicubic interpolation is handled below.
+  ## Limits indices to the bounds 1 and end.
+  if (strcmp(method, "linear"))
+    idx_rows = max (min (idx_rows, in_rows), 1);
+    idx_cols = max (min (idx_cols, in_cols), 1);
+  endif
 
   ## perform the actual resizing
   if (in_rows == out_rows && in_cols == out_cols)
     ## no resizing to do
-  elseif (strcmpi (method, "nearest"))
-    im = im(round(idx_rows), round(idx_cols), :);
-
+  elseif (strcmp (method, "nearest"))
+    im = im(round (idx_rows), round (idx_cols), :);
+  elseif (strcmp (method, "cubic"))
+    ## When rounding the output size up, it can happen that some interpolation points lap out at the right or the bottom of the input image.
+    ## In these cases we do not want to set the outlying points to 0 afterwards, but use the padded values. "symmetric" is also used in Matlab.
+    ## For an example, see tests below with scale 1/3.
+    im = imremap (im, idx_cols, idx_rows, method, "symmetric");
   else
-    idx_cols = repmat (idx_cols, [out_rows 1]);
-    idx_rows = repmat (idx_rows', [1 out_cols]);
     im = imremap (im, idx_cols, idx_rows, method);
   endif
 endfunction
@@ -265,15 +269,14 @@ endfunction
 %!         16   41   77  123  130  123  115  157  204  214  145   69   48   71   98
 %!         69   62   61   89  143  174  112  146  202  235  147   46   30   53   80
 %!        129   95   60   59  151  213   94  123  192  238  153   52   36   59   86];
-%! assert (max (max (abs (imresize (uint8 (in), 1.5, "bilinear") - uint8 (out)))) <= 1)
+%! assert (imresize (uint8 (in), 1.5, "bilinear"), uint8 (out), 1)
 %!
 ## Factor 0.5, gives an output of size 4 x 5 (without requiring to round the size)
 %! out = [ 98  152  126   57  112
 %!        129  124  171  145   83
 %!        116   95   80  116  164
 %!         61  143  146  146   53];
-%!
-%! assert (max (max (abs (imresize (uint8 (in), 0.5, "bilinear") - uint8 (out)))) <= 1)
+%! assert (imresize (uint8 (in), 0.5, "bilinear"), uint8 (out), 1)
 %!
 ## Factor 4/3, gives an output of size 10.6667 x 13.3333 rounded up to 11 x 14
 %! out = [116  185  199  143   80  117  164  177   76  103  133  135  110  106
@@ -287,20 +290,21 @@ endfunction
 %!         15   44   88  133  129  121  149  204  219  124   55   44   85   92
 %!         84   66   67  102  189  132  127  198  237  123   42   34   74   81
 %!        129   86   60   82  220  131  108  190  241  128   47   39   79   86];
-%! assert (max (max (abs (imresize (uint8 (in), 4/3, "bilinear") - uint8 (out)))) <= 1)
+%! assert (imresize (uint8 (in), 4/3, "bilinear"), uint8 (out), 1)
 %!
 ## Factor 1/3, gives an output of size 2.6667 x 3.3333 rounded up to 3 x 4
 %! out = [ 47   90   24   44
 %!        143   51   21   90
 %!         60   64   56   86];
-%! assert (max (max (abs (imresize (uint8 (in), 1/3, "bilinear") - uint8 (out)))) <= 1)
-
+%! assert (imresize (uint8 (in), 1/3, "bilinear"), uint8 (out), 1)
 
-## Comparison to the reference result from matlab with cubic interpolation
-## fails. This has to be fixed in imremap and here (remove the work around
-## above).
-%!xtest
-%!
+## Test bicubic interpolations against some reference results from matlab
+## without anti-aliasing with
+## imresize (uint8 (in), <FACTOR>, 'bicubic', 'Antialiasing', false)
+## The floating point error is less than 4e-13, but with integer Matlab rounds
+## wrong. So a difference of 1 is acceptable.
+## Factor 1.5, gives an output of size 12 x 15 (without requiring to round the size)
+%!test
 %! in = [116  227  153   69  146  194   59  130  139  106
 %!         2   47  137  249   90   75   16   24  158   44
 %!       155   68   46   84  166  156   69  204   32  152
@@ -322,8 +326,34 @@ endfunction
 %!          3   29   73  131  136  120  116  162  214  229  147   54   35   62   96
 %!         67   54   51   83  153  187  111  141  210  255  149   22   13   42   74
 %!        142   99   53   43  164  237   77  103  197  254  159   42   31   59   91];
-%! assert (imresize (uint8 (in), 1.5, "bicubic"), uint8 (out))
-
+%! assert (imresize (uint8 (in), 1.5, "bicubic"), uint8 (out), 1)
+%!
+## Factor 0.5, gives an output of size 4 x 5 (without requiring to round the size)
+%! out = [ 92   164   123    38   118
+%!        139   116   188   167    69
+%!        121    87    67   108   180
+%!         54   153   141   149    42];
+%! assert (imresize (uint8 (in), 0.5, "bicubic"), uint8 (out), 1)
+%!
+## Factor 4/3, gives an output of size 10.6667 x 13.3333 rounded up to 11 x 14
+%! out = [116   203   221   141    62   110   180   191    70   104   143   136   111   106
+%!         26    78   126   156   200   139   103   103    33    28    92   158    67    46
+%!         51    35    51   112   195   146   101    87    29    57   100   114    81    74
+%!        159   110    63    50    82   148   179   152    83   173   147    27   143   170
+%!         70   171   189   134   217   226   193   158   186   157    83    25   114   135
+%!         91   152   162   116   118   114   102    98   138    65    60   127    92    84
+%!         90   130   144   111    52    24    50   101    94    86   129   190   146   135
+%!         41   114   157   139    89    37    82   178   125   192   203   154   213   227
+%!          4    33    89   141   127   118   151   213   232   119    35    34    81    92
+%!         88    61    54    97   203   129   115   203   255   119    18    24    70    81
+%!        147    91    43    68   247   125    80   191   255   130    33    37    83    94];
+%! assert (imresize (uint8 (in), 4/3, "bicubic"), uint8 (out), 1)
+%!
+## Factor 1/3, gives an output of size 2.6667 x 3.3333 rounded up to 3 x 4
+%! out = [ 47    90    24    44
+%!        143    51    21    90
+%!         60    64    56    86];
+%! assert (imresize (uint8 (in), 1/3, "bicubic"), uint8 (out), 1)
 
 ## Reduce the size of a checkerboard, such that it turns into gray
 %!test
@@ -337,29 +367,59 @@ endfunction
 %!test
 %!
 %! in = ones (2) + 1i;
-%! out = imresize (in, 1.5);
-%! assert (out, ones (3) + 1i);
+%! out_nearest = imresize (in, 1.5, "nearest");
+%! assert (out_nearest, ones (3) + 1i);
+%! out_linear = imresize (in, 1.5, "linear");
+%! assert (out_linear, ones (3) + 1i);
+%! out_cubic = imresize (in, 1.5, "cubic");
+%! assert (out_cubic, ones (3) + 1i, 1e-14);
 
 ## check resizing of 1 pixel rgb images
 %!test
 %!
 %! in = cat (3, 10, 10, 10);
-%! out = imresize (in, [2, 2], "nearest");
-%! assert (size (out), [2, 2, 3]);
-%! out2 = imresize (in, [2, 2]);   # cubic interpolation (default)
-%! assert (size (out2), [2, 2, 3]);
+%! expected = 10 * ones (2, 2, 3);   # consistent with MATLAB behaviour
+%!
+%! out_nearest = imresize (in, [2, 2], "nearest");
+%! assert (out_nearest, expected);
+%! out_linear = imresize (in, [2, 2], "linear");
+%! assert (out_linear, expected);
+%! out_cubic = imresize (in, [2, 2], "cubic");
+%! assert (out_cubic, expected);
 
-## check resizing of line and row rgb images
+## check resizing of row and col rgb images
 %!test
 %!
-%! in = cat(3, [10,6],  [10,6],  [10,6]);
-%! out = imresize (in, [3, 3], "linear");
-%! plane_expected = [1; 1; 1] * [10, 8, 6];
-%! out_expected = cat (3, plane_expected, plane_expected, plane_expected);
-%! assert (out, out_expected);
+%! in_row = cat(3, [10,6],  [10,6],  [10,6]);
+%!
+%! plane_expected_nearest = [1; 1; 1] * [10, 10, 6, 6];
+%! out_expected = cat (3, plane_expected_nearest, plane_expected_nearest, plane_expected_nearest);
+%! out_nearest = imresize (in_row, [3, 4], "nearest");
+%! assert (out_nearest, out_expected);
+%!
+%! plane_expected_linear = [1; 1; 1] * [10, 8, 6];
+%! out_expected = cat (3, plane_expected_linear, plane_expected_linear, plane_expected_linear);
+%! out_linear = imresize (in_row, [3, 3], "linear");
+%! assert (out_linear, out_expected);
+%!
+%! plane_expected_cubic = [1; 1; 1] * [10.27777777777777, 8, 5.72222222222222];    # values from MATLAB
+%! out_expected = cat (3, plane_expected_cubic, plane_expected_cubic, plane_expected_cubic);
+%! out_cubic = imresize (in_row, [3, 3], "cubic");
+%! assert (out_cubic, out_expected, 1e-13);
 %!
-%! in2 = cat(3, [10;6],  [10;6],  [10;6]);
-%! out2 = imresize (in2, [3, 3], "linear");
-%! plane2_expected = [10; 8; 6] * [1, 1, 1];
-%! out2_expected = cat (3, plane2_expected, plane2_expected, plane2_expected);
-%! assert (out2, out2_expected);
+%! in_col = cat(3, [10;6],  [10;6],  [10;6]);
+%!
+%! plane_expected_nearest = [10; 10; 6; 6] * [1, 1, 1];
+%! out_expected = cat (3, plane_expected_nearest, plane_expected_nearest, plane_expected_nearest);
+%! out_nearest = imresize (in_col, [4, 3], "nearest");
+%! assert (out_nearest, out_expected);
+%!
+%! plane_expected_linear = [10; 8; 6] * [1, 1, 1];
+%! out_expected = cat (3, plane_expected_linear, plane_expected_linear, plane_expected_linear);
+%! out_linear = imresize (in_col, [3, 3], "linear");
+%! assert (out_linear, out_expected);
+%!
+%! plane_expected_cubic = [10.27777777777777; 8; 5.72222222222222] * [1, 1, 1];    # values from MATLAB
+%! out_expected = cat (3, plane_expected_cubic, plane_expected_cubic, plane_expected_cubic);
+%! out_cubic = imresize (in_col, [3, 3], "cubic");
+%! assert (out_cubic, out_expected, 1e-13);
