gauss sedel

function x = gauss_seidel(A, b, x0, k)
    D=diag(diag(A));
    L=tril(A)-D;
    U=triu(A)-D;
    x=x0;
    for i=1:k
        x=(D+L)^-1*(b-U*x);
    end
end