Here is the corrected code:
A = [2 0 1; -2 4 1; -1 -1 3];
b = [8 0 2]';
sz = size(A, 1);
disp([A b]);
% Forward elimination
for i = 1:sz-1
for j = i+1:sz
k = A(j, i) / A(i, i);
A(j, :) = A(j, :) - k * A(i, :);
b(j) = b(j) - k * b(i);
disp([A b]);
pause(3);
end
end
% Back Substitution:
x = zeros(sz, 1);
for i = sz:-1:1
x(i) = (b(i) - A(i, i+1:end) * x(i+1:end)) / A(i, i);
end
disp('FINAL Solution: ');
disp(x);