% Testing the smoothing properties of the Jacobi iteration
% 1D Laplace
global V DA
global y_ex

n=input('Enter matrix size: ');
h=1/n; % mesh-size
x=[h:h:1-h];
I=[1:1:n-1];
% Creat 1D Laplacian
e = ones(n,1);
A = spdiags([-e 2*e -e], -1:1, n-1, n-1);
[V,DA]=eig(full(A));
[V,DA]=sorteer(V,DA);

y0  = zeros(n-1,1);
y_ex= x.^2.*sin(x);y_ex=y_ex';
b   = A*y_ex;
eps_jac = 1e-3;
flag=0;

disp('Run Jacobi iteration method...')
[y,r,it_Jac] = Jacobi(A,b,y0,eps_jac,flag);

disp(['Converged to ' num2str(eps_jac) ' for ' int2str(it_Jac) ' iter.'])
disp('Rerun and trace the error components...') ,wait

flag=1;
[y,r,it_Jac] = Jacobi(A,b,y0,eps_jac,flag);

return


plot(4*(sin(I*pi/(2*n))).^2,'or'),hold
plot(diag(DA)) 

d=diag(A);
di = 1./d;
D=spdiags(d,0,n-1,n-1);
R=D-A;
B=inv(D)*R;
EB=eig(full(B));
rho=max(abs(EB));