How can we analytically derive the flow-lines of a normal permanent bar-magnet?
Physics context & own approach:
In classical electromagnetics we have the legendary Maxwell's Equations:
$$\begin{align}\nabla \cdot {\bf E} &= \frac{\rho}{\epsilon_0}\\\nabla \cdot {\bf B} &=0\\\nabla \times {\bf E} &= -\frac{\partial \bf B}{\partial t}\\\nabla\times {\bf B} &= \mu_0 \left( {\bf J}+\epsilon_0\frac{\partial {\bf E}}{\partial t} \right)\end{align}$$
The situation in a normal permanent magnet far away from any electric field, we should have:
$$\begin{align}{\bf E} &= 0\\\bf J &= 0\\{\bf B} &= \cases{c{\bf \hat y} \hspace{1.0cm}\text{(inside magnet)} \\{\bf B}(x,y) \text{ (outside)}}\end{align}$$ (c is constant : homogeneous magnet with constant microscopic dipole moment)
Can we use this to derive analytic expression for the flow-lines (or vector field $\bf B$) around a permanent magnet?
Pictures show color mapped representation of B field around bar and horseshoe magnets acquired from a numerical optimization.
Here is what I suspect is a standard approach to such problems:
We have
$\nabla \times \mathbf B = 0; \tag 1$
then
$\nabla \times (\nabla \times \mathbf B) = 0; \tag 2$
there is a standard identity from vector calculus which asserts that
$\nabla \times (\nabla \times \mathbf B) = \nabla(\nabla \cdot \mathbf B) - \nabla^2 \mathbf B; \tag3$
using this in concert with the Maxwell equation
$\nabla \cdot \mathbf B = 0 \tag 4$
yields
$\nabla \times (\nabla \times \mathbf B) = - \nabla^2 \mathbf B; \tag5$
(2) thus becomes
$\nabla^2 \mathbf B = 0, \tag 6$
where it is understood that $\nabla^2$ is applied to $\mathbf B$ component-wise; thus we have
$\nabla^2 \mathbf B_x = \nabla^2 \mathbf B_y = \nabla^2 \mathbf B_z = 0; \tag 7$
we may impose boundary conditions on the components of by by assuming
$\mathbf B = c \hat{\mathbf y} \tag 8$
in the iron bar, as is given in the text of the problem. We also assume that
$\mathbf B \to 0 \; \text{at infinity}; \tag 9$
since the components of $\mathbf B$ satisfy Laplace's equation, this should be sufficient information to assure the existence of a unique solution for $\mathbf B$; but I strongly suspect it will not be very nice, analytically speaking.
To conclude, there is probably no "nice" formula for the field lines, but they can of course be plotted via numerical methods once we have found the vector field $\mathbf B$.