clear, clc;

U1 = zeros(3,3);

for i=1:3
    for j=1:3
        if (j >= i)
            U1(i,j) = 1;
        end
    end
end

b = [1; 2; 3];

x = backward_substitution(U1,b);

x

U1\b

norm(b-U1*x,2)

