Update:
After a long struggle,
It seems this might just be solving the Laplace equation (zero source coulomb gauge!) in cylindrical coordinates, computing the components of B from $\phi$, and then Taylor expanding.
I am revisiting problem 5.4 in Jackson's electrodynamics
A magnetic induction B in a current-free region in a uniform medium is cylindrically symmetric with components $B_z(\rho,z)$ and $B_{\rho}(\rho,z)$ and with a known $B_z(0,z)$ on the axis of symmetry. The magnitude of the axial field varies slowly in z.
Show that near the axis the axial and radial components of magnetic induction are approximately
$B_z (\rho,z) \approx B_z(0,z) - (\frac{\rho^2 }{4} ) [ \frac{\partial^2 B_z(0,z)}{ \partial^2 z }] + . . .$
and
$B_{\rho}( \rho,z) \approx -(\rho/2) [ \frac{\partial B_z(0,z)}{ \partial z}] + (\frac{\rho^3}{16}) [ \frac{\partial^3 B_z(0,z)}{ \partial z^3}] + . . ..$
