function x = SOR(A,b,x0,w)
  
L = -tril(A,-1);
U = -triu(A,1);
D = diag(diag(A));

Nmax = 500;

xe = [-1;1;-1;1];
x = x0;
for j=1:Nmax
  
  x = (D-w*L)\((w*U+(1-w)*D)*x+w*b);
  
  if max(abs(x-xe)) < 10^(-4)
    break
  end
    
end
iters = j