nx = 63;
[~,~,A]= laplacian([nx,nx]);


N = size(A,1);
x = zeros(N,1);
b = ones(N,1) / (N+1)^2;
norm_b = norm(b);

max_iter = 10000;
res = NaN(1,max_iter);
tic
for j=1:max_iter
	x = gauss_seidel(N, A,x,b);
	norm_res = norm(b -A*x);
	res(j) = norm_res;
	if ( norm_res/norm_b < 1.0e-10 )
		break;
		end
end
T = toc;
disp([num2str(j),' iterations,   elapsed time: ',num2str(T),' sec.']);


semilogy( 1:j, res(1:j), 'bx-', 'LineWidth',2 )
axis tight
xlabel('# iterations')
ylabel('|| b -A x^j ||')
