The issue with your approach is that you are trying to use matrix addition between two arrays with different shapes. The first loop results in a 3D matrix of size (2 * range + 1)^3, whereas the second approach results in a 1D vector of size (2 * range + 1). You cannot add these two arrays together because they have different shapes.
To solve this problem, you can try to reshape the 1D vector into a 3D matrix of the same size as the output of the first loop before adding it to the canvas matrix. The following code should work:
[x, y, z] = ndgrid(-range:1:range, -range:1:range, -range:1:range);
result = const * exp(-sum(((delta * [x(:), y(:), z(:)] * RotationMatix) ./ [a, a, b]).^2));
%reshape so the result can be added to the canvas
result = reshape(result, size(x));
canvas(index(1)-range:index(1)+range, index(2)-range:index(2)+range, index(3)-range:index(3)+range) += result;
