--- a/scripts/general/cumtrapz.m
+++ b/scripts/general/cumtrapz.m
@@ -50,53 +50,61 @@
     print_usage ();
   endif
 
-  nd = ndims (x);
-  sz = size (x);
+  have_xy = have_dim = false;
 
-  have_x = false;
-  have_dim = false;
   if (nargin == 3)
-    have_x = true;
+    have_xy = true;
     have_dim = true;
   elseif (nargin == 2)
     if (! size_equal (x, y) && isscalar (y))
       dim = y;
       have_dim = true;
     else
-      have_x = true;
+      have_xy = true;
     endif
   endif
 
+  if (have_xy)
+    nd = ndims (y);
+    sz = size (y);
+  else
+    nd = ndims (x);
+    sz = size (x);
+  endif
+
   if (! have_dim)
     ## Find the first non-singleton dimension.
-    dim = find (sz > 1, 1);
-    if (isempty (dim))
-      dim = 1;
-    endif
+    (dim = find (sz > 1, 1)) || (dim = 1);
   else
-    dim = floor (dim);
-    if (! (isscalar (dim) && 1 <= dim && dim <= nd))
-      error ("cumtrapz: invalid dimension DIM");
+    if (!(isscalar (dim) && dim == fix (dim))
+        || !(1 <= dim && dim <= nd))
+      error ("trapz: DIM must be an integer and a valid dimension");
     endif
   endif
 
   n = sz(dim);
-  idx1 = cell ();
-  for i = 1:nd
-    idx1{i} = 1:sz(i);
-  endfor
-  idx2 = idx1;
+  idx1 = idx2 = repmat ({:}, [nd, 1]);
   idx1{dim} = 2 : n;
   idx2{dim} = 1 : (n - 1);
 
-  if (! have_x)
+  if (! have_xy)
     z = 0.5 * cumsum (x(idx1{:}) + x(idx2{:}), dim);
   else
-    if (! size_equal (x, y))
-      error ("cumtrapz: X and Y must have the same shape");
+    if (isvector (x) && !isvector (y))
+      if (length (x) != sz(dim))
+        error ("cumtrapz: length of X and length of Y along DIM must match");
+      endif
+      ## Reshape vector to point along dimension DIM
+      shape = ones (nd, 1);
+      shape(dim) = sz(dim);
+      x = reshape (x, shape);
+      z = 0.5 * cumsum (bsxfun (@times, diff (x), y(idx1{:}) + y(idx2{:})), dim);
+    else
+      if (! size_equal (x, y))
+        error ("cumtrapz: X and Y must have same shape");
+      endif
+      z = 0.5 * cumsum (diff (x, 1, dim) .* (y(idx1{:}) + y(idx2{:})), dim);
     endif
-    z = 0.5 * cumsum ((x(idx1{:}) - x(idx2{:})) .*
-                      (y(idx1{:}) + y(idx2{:})), dim);
   endif
 
   sz(dim) = 1;
@@ -104,13 +112,23 @@
 
 endfunction
 
+
 %!shared x1,x2,y
 %! x1 = [0,0,0;2,2,2];
 %! x2 = [0,2,4;0,2,4];
 %! y = [1,2,3;4,5,6];
-%!assert (cumtrapz(y),[0,0,0;2.5,3.5,4.5])
-%!assert (cumtrapz(x1,y),[0,0,0;5,7,9])
-%!assert (cumtrapz(y,1),[0,0,0;2.5,3.5,4.5])
-%!assert (cumtrapz(x1,y,1),[0,0,0;5,7,9])
-%!assert (cumtrapz(y,2),[0,1.5,4;0,4.5,10])
-%!assert (cumtrapz(x2,y,2),[0,3,8;0,9,20])
+%!assert (cumtrapz(y), [0,0,0;2.5,3.5,4.5])
+%!assert (cumtrapz(x1,y), [0,0,0;5,7,9])
+%!assert (cumtrapz(y,1), [0,0,0;2.5,3.5,4.5])
+%!assert (cumtrapz(x1,y,1), [0,0,0;5,7,9])
+%!assert (cumtrapz(y,2), [0,1.5,4;0,4.5,10])
+%!assert (cumtrapz(x2,y,2), [0,3,8;0,9,20])
+%% Test ND-array implementation
+%!shared x1,x2,y
+%! x1 = 1:3;
+%! x2 = reshape ([0,2,4;0,2,4], [1 2 3]);
+%! y = reshape ([1,2,3;4,5,6], [1 2 3]);
+%!assert (cumtrapz(y,3), reshape([0,1.5,4;0,4.5,10],[1 2 3]))
+%!assert (cumtrapz(x1,y,3), reshape([0,1.5,4;0,4.5,10],[1 2 3]))
+%!assert (cumtrapz(x2,y,3), reshape([0,3,8;0,9,20],[1 2 3]))
+