% Simple numerical indefinite integral
function result = integral(integrand, dx)

result = integrand; % set up result with same dimensions as integrand

% Shift so that the first element of the integral will be zero.
i2 = circshift(vec(integrand),1);
   % circshift requires that the argument be a column vector.
i2(1) = 0;
result(:) = cumsum(i2)*dx; %result is filled with the cumulative sum, but
   % its dimensions remain the same as the original integrand array.