"""
    laplacian(U::Matrix{Float64}, N::Int, h::Float64)

    Computes the discrete 2D Laplacian of a matrix `U` using a 5-point stencil
    and circular boundary conditions.

    # Arguments
    - `U::Matrix{Float64}`: The input 2D matrix representing the field or image.
    - `N::Int`: Integer
    - `h::Float64`: The spatial step size or grid spacing between points in the discretization.

    # Returns
    - `Vector{Float64}`: A flattened (vectorized) representation of the approximated Laplacian values for each element in `U`. The boundary conditions are handled circularly.
"""
function laplacian(U::Matrix{Float64}, N::Int, h::Float64)
    Δ = zeros(N, N)
    for i in 2:N-1, j in 2:N-1
        Δ[i, j] = (U[i+1, j] + U[i-1, j] + U[i, j+1] + U[i, j-1] - 4 * U[i, j]) / h^2
    end
    return vec(Δ)
end