h = 1/64;
nx = 1/h -1;
[~,~,A_h]= laplacian([nx,nx]);
n = size(A_h,1);

H = 2*h;
NX = 1/H -1;
[~,~,A_H]= laplacian([NX,NX]);
N = size(A_H,1);


x_h = zeros(n,1);
b_h = ones(n,1) / (n+1)^2;
norm_b = norm(b_h);


nu_1 = 3;
nu_2 = 3;

max_iter = 10000;
res = NaN(1,max_iter);
tic
for j=1:max_iter
	x_h = twogrid(n, A_h,x_h,b_h,   nu_1,nu_2, N,A_H,nx,NX);
	norm_res = norm(b_h -A_h*x_h);
	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 ||')
