function laplacian(u, N, h)
    U = reshape(u, N, N)
    padded = circshift(U, (-1, 0)) .+ circshift(U, (1, 0)) .+
             circshift(U, (0, -1)) .+ circshift(U, (0, 1)) .- 4 .* U
    return vec(padded) ./ h^2
end