function x = jacobi(A,b,x0,tol,MAX_ITER)

%%%%%
% This is a Jacobi method to solve Ax = b using the max norm for tolerance
%%%%%

N = length(b);
c = zeros(N,1);
T = zeros(N);
x = zeros(N,1);
strformat = ['%d'];
for i=1:N,
    c(i) = b(i)/A(i,i);
    T(i,:) = -A(i,:)/A(i,i);
    T(i,i) = 0;
    strformat = [strformat,'  %1.4f  '];
end
strformat = [strformat,' \n'];

k = 0;
fprintf('---------------------------------\n');
fprintf('  k              x           \n');
fprintf('---------------------------------\n');
while( k < MAX_ITER )
    fprintf(1,strformat,k,x0);
    x = T*x0' + c;
    test = max(abs(x0-x'));
    x0 = x';
    k = k+1;
    if(test < tol)
        fprintf(1,strformat,k,x0);
        return;
    end    
end