function [X,iter,err] = JACOBI(A,B,epsilon,Niter,tol); 
% [X,iter,err] = JACOBI(A,B,epsilon,Niter,tol); 
[N,N] = size(A); DA = diag(diag(A)); EA = DA - A; 
IDA   = DA^(-1); AA = IDA*EA; 
X0    = 2*tol*(rand(N,1)-1/2); X = X0;  
err   = 2*epsilon; iter = 0;   
while (err > epsilon) && (iter < Niter)
   Xold = X; X = AA*X + IDA*B; 
   err  = norm(X-Xold); iter = iter + 1;
end
end