function [y_next,r,k] = Jacobi(A,b,y0,eps_jac,flag)

global V DA
global y_ex

figure(1),clf,hold on

n=size(A,1);
y_prev=y0;
k=1;
r(k)=norm(b);

d=diag(A);
di = 1./d;
D=spdiags(d,0,n,n);
R=D-A;
while (r(k) > eps_jac),
      y_next=di.*(R*y_prev+b);
      k=k+1;
      r(k)=norm(b-A*y_next);
      if flag==1,
         err=y_ex-y_next;
         erc=err'*V;
         plot(erc,'dr'),wait
      end
      y_prev=y_next;
end
k=k-1;
