function [A,b] = poisson2(N)

e = ones(N,1);
e_small = ones(N-1,1);

B = 4*diag(e) - diag(e_small,-1) - diag(e_small,1);
I = -diag(e);

A = kron(diag(e),B) + kron(diag(e_small,1),I) + kron(diag(e_small,-1),I);

h = 1 / (N + 1);

Y = h * meshgrid(1 : N);
X = Y';

B = 5 * pi^2 * sin(2 * pi * X) .* sin(pi * Y);
b = -h^2*reshape(B, N^2, 1);