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

function[x,s]=jacobi(A,B,X0,eps)
n=length(B);
D=diag(A);
R=A-diag(D);
x=X0; s=0; xreal=ones(n,1);
forwarderr=x-xreal;
while norm(forwarderr,inf)>eps
    x=(B-R*x)./D; 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]=jacobi(A,B,X0,eps)

x =

   0.999999779847827
   0.999999578485665
   0.999999409667645
   0.999999290421916
   0.999999227328151
   0.999999227328151
   0.999999290421916
   0.999999409667645
   0.999999578485665
   0.999999779847827


step =

    32