function    [x,X]=jacobi_2(A,b,x0,N,TOL);
%function    [x,X]=jacobi_2(A,b,x0,N,TOL);
%%%%%%%%%%%%%%%%%%% Jacobi Iterative Algorithm %%%%%%%%%%%%%%%%%%%%%%
%%     
%%     Solve Ax=b by Jacobi iteration with initial approximation x0
%%   INPUTS:
%%     A -  the coefficent matrix of the system
%%     b -  the "Right Hand Side'' of the equation
%%     N - maximum number of iterations permitted
%%     TOL - tolerance
%%     
%%   OUTPUT:
%%     x - approximate solution or a message to say max number of
%%                                         iterations is exceeded.
%%     X - all the iterations (each as a row)

% Step 1
k=1;
D=diag(A);DD=diag(ones(D)./D);
U=triu(A-diag(diag(A)));
L=tril(A-diag(diag(A)));
flag=0;
% Step 2
while  (k<=N & flag==0);
    x=-DD*(L+U)*x0+DD*b;
   if norm(x-x0)< TOL;
      flag=1;
   end;
   X=[X;x'];
   x0=x;k=k+1;
end;
if flag==0;
   message = 'maximum number of iterations exceeded'
end;
