function x=matrixSOR(A,b,n,w)
D=diag(diag(A));
L=-tril(A,-1);
U=-triu(A,1);
BSOR=inv(D-w*L)*[(1-w)*D+w*U];
f=w*inv(D-w*L)*b;
x=zeros(n,1);
rx=BSOR*x+f;
while norm(rx-x)>=1.0e-6
    x=rx;
    rx=BSOR*x+f;
end
x=rx;