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

function[x s]=GaussSeidel(A,B,X0,eps)
% 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(L+D)*(-U*x+B); forwarderr=x-xreal;
    s=s+1;
end
___________________________________________________________________________
format long
n=10; eps=10^(-6);
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]=GaussSeidel(A,B,X0,eps)

x =

   0.999999415556504
   0.999999273676315
   0.999999342118390
   0.999999487297440
   0.999999639349225
   0.999999767354707
   0.999999862406187
   0.999999926597182
   0.999999966352989
   0.999999988784330


step =

    18