function laplacian(U::Matrix{Float64}, N::Int, h::Float64)
    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