function [A,b] = twodlaplace(N)
% 2D Laplace auf dem Einheitsquadrat
% mit rechter Seite x^3*(x-1)*y*(y-1)

h = 1/N;
A = gallery('poisson',N-1);
x = (h:h:(N-1)*h)';

% ! Funktioniert nur so einfach weil die
% ! rechte Seite separierbar ist.
b = h^2*reshape((x.^3.*(x-1))*(x.*(x-1))',(N-1)^2,1);