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))