function [x]=jacobi(A,b,x0,N)
n=length(b);
x=x0;
for k=1:N,
    xs=zeros(n,1);
    for i=1:n,
        s=b(i);
        for j=1:n
            if i~=j,
                s=s-A(i,j)*x(j);
            end;
        end;
        xs(i)=s/A(i,i);
    end;
    x=xs;
end;
