function [x,iter] = jacobi(A,b,tol,maxit)
% uses Jacobi iteration to solve A*x=B
% A      square  n x n
% b      n x 1
% tol    tolerance
% maxit  maximal number of iterations
% iter   actual number of iterations required for convergence

[n,m]=size(A)
x = zeros(size(b));      %initial guess
i=1;
iter = maxit;
for i = 1: maxit
    for j=1:n
        y(j)=(b(j)-A(j,1:j-1)*x(1:j-1)-A(j,j+1:n)*x(j+1:n))/A(j,j);
    end
    if max(abs(A*y'-b)) < tol
        iter = i;
	x=y';
        break
    end
    x=y';
end
