clear;
clf;

N = 16;

A = (N+1)^2*(2 * eye(N) - diag(ones(N-1,1), -1)  - diag(ones(N-1,1), 1));
D = diag(diag(A));
E = tril(A);  


% Let's generate a random right hand side

b = 2*ones(N,1);
b(N /4+1:3*N/4) = -1;
x_true = A\b;


%x_true = sin( 7*pi * (1:N)'/(N+1)) + sin( 2*pi * (1:N)'/(N+1));
%b = A * x_true;

% True solution
plot(x_true)

xJ = zeros(N,1);
xGS = zeros(N,1);


subplot(1,3,1)
plot((0:N+1)'/(N+1),[0;x_true - xJ;0]);
title('Error in Jacobi');
eJ(1) = norm(x_true);
rJ(1) = norm(b);
h = gca();
ylimJ = h.YLim;


subplot(1,3,2)
plot((0:N+1)'/(N+1),[0;x_true - xGS;0]);
title('Error in Gauss-Seidel');
eGS(1) = norm(x_true);
rGS(1) = norm(b);

pause
%wj = 2/3;
%wgs = 1.5;
for i = 2:100
%    pause
    xJ = xJ + D \ (b - A*xJ);
    xGS = xGS + E \ (b - A*xGS);
    
  
    subplot(1,3,1)
    plot((0:N+1)'/(N+1),[0;x_true - xJ;0]);
    ylim(ylimJ);
    title('Error in Jacobi');
    
    eJ(i) = norm(x_true-xJ);
    
    rJ(i) = norm(b - A*xJ);
    
    subplot(1,3,2)
    plot((0:N+1)'/(N+1),[0;x_true - xGS;0]);
    ylim(ylimJ);
    title('Error in Gauss-Seidel');
    
    eGS(i) = norm(x_true-xGS);
    
    rGS(i) = norm(b-A*xGS);
    subplot(1,3,3)
    %plot([rJ', rGS'])
    semilogy([rJ', rGS'])
    legend('Jacobi', 'G-S')
    title('Residuals')
    pause(.1)
end
