% SOR Method solves AX=B.
% Input: A - nxn coefficient matrix, B - nx1 vector 
%        X0 - nx1 initial vector, eps - tolerance (forward), w - weight 

function[x s]=SOR(A,B,X0,eps,w)
% Part: A=L+D+U
n=length(B);
L=zeros(n,n);
U=zeros(n,n);
for j=1:n-1
for k=j+1:n
L(k,j)=A(k,j);
end
end
D1=diag(A);
D=diag(D1);
U=A-D-L;
% end
% Iteration part
x=X0; s=0; xreal=ones(n,1);
forwarderr=x-xreal;
while norm(forwarderr,inf)>eps
    x=inv(w*L+D)*(-w*U*x+w*B+(1-w)*D*x); forwarderr=x-xreal;
    s=s+1;
end
___________________________________________________________________________
format long
n=10; eps=10^(-6); w=1.2;
e=ones(n,1); X0=zeros(n,1);
A=spdiags([-e 3*e -e],-1:1,n,n); % Constructing A matrix.
B=ones(n,1); B(1)=2; B(n)=2;     % Adjusting two entries of B.
[x step]=SOR(A,B,X0,eps,w)

x =

   1.000000329753529
   1.000000178833364
   1.000000054803257
   1.000000007116585
   0.999999996167126
   0.999999996096381
   0.999999997674279
   0.999999998809740
   0.999999999697792
   0.999999999299867


step =

    13