function [x,J,c] = jacobi(A,b,n,z)
%
%   x = jacobi(A,b,n,z)
%    
%      Jacobi iteration on system A*x = b with printing
%      n -- number of iterations
%      z -- initial vector  (default 0)
%
%      x -- final iterate
%      J -- Jacobi matrix
%      c -- Jacobi vector
%

if nargin <=3, z=0*b; end

D = diag(diag(A));
J = D\(D - A);
c = D\b;
x=z;
for k = 1:n
   x = J*x + c;
   fprintf(1,'%3d     ',k)
   fprintf(1,'%5.4f     ',x')
   fprintf(1,'\n')
end