For loop vs. block solution does not nearly same result

I am solving trying to solve the following sparse matrix system, initially as follows. gg_n0 is a vector of length I*J*L.

Q = (speye(I*J*L) - AT*Delta);
gg_n1 = Q\gg_n0; % original

Because the block diagonal Q matrix is quite large, I try to loop it which is much faster in this case.

for l=1:L 
    gg_n1_alt((l-1)*I*J+1:l*I*J)= Q((l-1)*I*J+1: l*I*J, (l-1)*I*J+1: l*I*J)\gg_n0((l-1)*I*J+1:l*I*J);
end

The problem is, as the original gg_n0 reshaped to it's marginal distributions in I and J (it is a probability distribution P(a,b,z) where a is discretized in I, b in J and z in L) sums to 1, the marginal distributions of gg_n1 should also sum to 1. In the first case of the simple mldivide this holds. However when doing the for loop, it does not.

Plotting the resulting gg_n1 by mldivide and no for loop yields

And with for loop (gg_n1_alt)

As such they are quite a like but still differ enough that the second is unusable. I suspect it may have to do with the fact that the numbers are so small?