function [out,Tw,cw] = sor(L,U,D,b,x0,n,w)
x = x0;
Tw = (D-w*L)\((1-w)*D+w*U);
cw = w*((D-w*L)\b);

for i=2:n
    x(:,i) = Tw*x(:,i-1) + cw;
end
out = x;
return;