[g,b] = sqgrid(n);
A=lapdir(g,b);
nn=sum(g(:)>0);
f=zeros(nn,1);
f(1:b)=sin((1:b)*4*pi/b);
u=A\f;
for i=1:n
    for j=1:n
        U(i,j)=u(g(i,j));
    end
end
mesh(U);