function [x_h] = twogrid(n, A_h,x_h,b_h,  nu_1,nu_2, N,A_H,nx,NX)

%  pre-smoothing iteration
for j=1:nu_1
	x_h = gauss_seidel(n, A_h,x_h,b_h);
	end

res_h = b_h -A_h*x_h;

%  restriction of residual to coarse grid
res_h = reshape(res_h,nx,nx)';
res_H =   res_h( 2:2:nx-1, 2:2:nx-1 ) ...
	+0.5* res_h( 1:2:nx-2, 2:2:nx-1 ) ...
	+0.5* res_h( 3:2:nx,   2:2:nx-1 ) ...
	+0.5* res_h( 2:2:nx-1, 1:2:nx-2 ) ...
	+0.5* res_h( 2:2:nx-1, 3:2:nx   ) ...
	+0.5* res_h( 1:2:nx-2, 1:2:nx-2 ) ...
	+0.5* res_h( 3:2:nx  , 3:2:nx   );
res_H = reshape(res_H' / 4, N,1);



%  %  exact solve on coarse grid
%  c_H = A_H \ res_H;

%  inexact base solver: 50 steps of iteration
c_H = zeros(N,1);
for j=1:50
	c_H = gauss_seidel(N, A_H,c_H,res_H);
	end


%  prolongation of coarse grid correction
c_H = reshape( c_H, NX,NX)';
c_h = zeros(nx);
c_h( 2:2:nx-1, 2:2:nx-1 ) = c_H;
c_h( 3:2:nx,   2:2:nx-1 ) = 0.5* c_H;
c_h( 2:2:nx-1, 3:2:nx   ) = 0.5* c_H;
c_h( 3:2:nx  , 3:2:nx   ) = 0.5* c_H;
c_h( 1:2:nx-2, 2:2:nx-1 ) = c_h( 1:2:nx-2, 2:2:nx-1 ) +0.5* c_H;
c_h( 2:2:nx-1, 1:2:nx-2 ) = c_h( 2:2:nx-1, 1:2:nx-2 ) +0.5* c_H;
c_h( 1:2:nx-2, 1:2:nx-2 ) = c_h( 1:2:nx-2, 1:2:nx-2 ) +0.5* c_H;
c_h = reshape(c_h',n,1);


x_h = x_h +c_h;


%  post-smoothing iteration
for j=1:nu_2
	x_h = gauss_seidel(n, A_h,x_h,b_h);
	end
