    N = 256;
    h = 1/N;
    
    % N-1 points per row and per per column
        
    e = ones(N-1,1);
    I = eye(N-1);
    K = spdiags([-e 2*e -e], [-1 0 1], N-1,N-1);
    L = (kron(I,K) + kron(K,I))/(h^2);
    
    % Full matrix
    % fullL = full(L);
    
    % Solve the Poisson equation, sparse direct solver
    f = ones(N-1,N-1);
    f = reshape(f,(N-1)^2,1);
    u = L \ f;
    u = reshape(u,N-1,N-1);
    
    figure(1); clf; 
    imagesc(u);