I've found this code very helpful. However, there is a minor issue. The output in the case of scalar arguments d1 and d2 is a column vector, while in the case of arguments of higher dimensionality, the output is transposed (it is a matrix consisting of raw vectors, corresponding to each element of the argument). If I'm not clear, just take a look how the output looks like when d1 and d2 are scalars, and when they are vectors 1 by 2.

I corrected the code by adding in the end of the function the following line:

if NDim == 0, y = y'; end

It changes the output of the scalar arguments case a raw vector, consistent with argument of higher dimensionality and with the original linspace function.

In 2D you can use these simple commands, that are much faster:
A = randi(10,1000,1);
B = randi(10,1000,1)+20;
N = 1500;
dx = (B-A)/(N-1);
AB = repmat(dx,1,N);
AB(:,1) = A;
AB = cumsum(AB,2);

I've found this code very helpful. However, there is a minor issue. The output in the case of scalar arguments d1 and d2 is a column vector, while in the case of arguments of higher dimensionality, the output is transposed (it is a matrix consisting of raw vectors, corresponding to each element of the argument). If I'm not clear, just take a look how the output looks like when d1 and d2 are scalars, and when they are vectors 1 by 2.
I corrected the code by adding in the end of the function the following line:
if NDim == 0, y = y'; end
It changes the output of the scalar arguments case a raw vector, consistent with argument of higher dimensionality and with the original linspace function.

Comment only