Some points to consider are:
-
Is the flux \mathbf{F} solenoidal (i.e., does \nabla\cdot\mathbf{F} = 0)? A vector field must be divergence-free to have a streamfunction representation. You would only expect a solenoidal heat flux in the absence of any distributed heat sources.
-
Solving the Poisson equation with pure Neumann boundary conditions is ill-posed, which can lead to numerical issues if not done carefully, as discussed, e.g., in this thread.