function [x, iter, time] = jac(A,b, x0, tol)

% jacobi to solve linear system A*x=b
% in:	A	matrix of linear system
%	b	right hand side of system
%	x0	initial guess
%	tol	tolerance for convergence
% out:	x	approximative solution 
%	iter	number of iterations needed
% 	time	cpu time needed
% auth: fab@math.ethz.ch

tic;		% start measuring time
D=diag(diag(A)); L=tril(A)-D; R=triu(A)-D;	% decompose matrix
S=-(L+R);
x=x0;	% start with initial guess
maxIter=2e3; iter=0;
while (iter<maxIter)
  iter=iter+1;
  xold=x;
  x=D\(S*x+b);	% Jacobi iteration without inverse
  if (norm(x-xold)/norm(x)<tol) break; end	% converged!
end
time=toc;	% stop measuring time
if (iter==maxIter) error(['no convergence after ',int2str(iter),' iterations']); end;