sor

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