
N = 2^7-1;
w = 1;  % w = 0 gives standard model problem  
h = 1/(N+1);
A = helmA(N, w);

[x,y] = meshgrid( (1:N)*h, (1:N)*h);

%To test, choose a solution and compute the associated 
%right hand side. 
uex = x .* (1-x) .* y .* (1-y);
uex=uex(:);
b = A*uex; 

surf(x,y,reshape(uex,N,N))

